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Many  current  fighter  aircraft  utilize  thin  delta  wings  with  forebody  strakes,  which 
produce  strong  vortical  flows  at  moderate  to  high  angles  of  attack.  The  resuiting  lift  incre¬ 
ment  is  highly  nonlinear  and  cannot  be  accurately  predicted  by  current  design  methods. 

The  objective  of  this  thesis  is  to  calculate  the  flow  over  a  cranked  delta  wing  using  the 
thin-layer  Navier-Stokes  equations.  Emphasis  is  placed  on  determining* the  effects  of  the 
strake  vortex  on  the  wing  vortex  and  the  ability  of  the  code  to  reproduce  the  secondary 
vortex  structure. 

The  algorithm  used  in  this  study  was  the  ARC3D  code,  written  by  Dr.  Thomas 
Pulliam  of  the  NASA  Ames  Research  Center.  It  has  been  extensively  modified  by  Dr. 
Philip  Webster  of  the  Flight  Dynamics  Laboratory.  The  calculations  were  done  using  the 
CDC  Cyber  and  Cray  XMP-12  computers  at  Wright-Patterson  AFB,  Ohio. 
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% 


Tabic  of  Cpntcms 


i 


i 


Preface . 

List  of  Figures . 

List  of  Symbols . 

Abstract . 

I.  Introduction . 

II.  Analysis . 

Governing  Equations . 

Coordinate  Transformation  .  .  .  . 
Thin-Layer  Approximation  .  .  .  . 
Boundary  Conditions . 

III.  Numerical  Procedure . 

Implicit  Time  Marching  Algorithm 

Linearization . 

Approximatate  Factorization.  .  .  . 

Diagonalization . 

Artificial  Dissipation . 

Solution  Procedure . 

SSD  Implementation . 

IV.  Results  and  Discussion . 

Delta  Wing . 

Cranked  Delta  Wing . 


Page 

11 

v 

vi 
viii 

1 

8 

8 

10 

12 
14 
16 
16 

17 

18 
19 
21 
22 

24 

25 
25 
40 
52 
54 


4 

5 

-1 


4 

4 

J 


1 


i 


V.  Conclusions  and  Recommendation . 

Appendix  A  :  Coordinate  Transformation  of  the  Governing  Equations 
Appendix  B  :  Flux  Jacobian  Matrices . 


59 


I 


Appendix  C  :  Flux  Jacobian  Eigenvalue  and  Eigenvector  Matrices .  62 

Appendix  D  :  Incore  Versus  Solid  State  Disk  Comparison .  66 

Appendix  E  :  Local  Application .  67 

Appendix  F  :  Algebraic  Grid  Generator .  70 

Bibliography .  84 

VITA .  88 


C 


C 


c 


ufet r-VjukV  4.  i»?«w  ■ . >•-..  •  ' 


Figure 

1  Vortex  Flow  Over  a  Delta  Wing 


Possible  Vortex  Flows  for  Sharp  Delta  Wings. 
75°  Delta  Wing  Geometry . 


Streamwise  Grid  Flane  for  75°  Delta  Wing 


Comparison  of  Pitot  Pressures  with  Experiment. 
Comparison  with  Full  Navier-Stokes  Solution  . 


Delta  Wing  Pitot  Pressures  and  Crossplane  Velocities 


8  Delta  Wing  Surface  Pressures 


80°/69°  Cranked  Delta  Wing  Geometry. 


10  Streamwise  Grid  Plane  for  80°/69°  Cranked  Delta  Wing  . 

1 1  Cranked  Wing  Pitot  Pressures  and  Crossplane  Velocities 


12  Comparison  of  Pressures  Upstream  and  Downstream  of  the  Crank. 

13  Cranked  Wing  Surface  Pressure  Comparison  With  Experiment.  .  . 


v 


a  speed  of  sound  (m/sec) 

Cp  specific  heat  at  constant  pressure,  or  coefficient  of  pressure 

E  in  viscid  flux  vector  in  <;  direction 

Ev  viscous  flux  vector  in  the  \  direction 

e  total  energy 

e,  specific  internal  energy 

F  inviscid  flux  vector  in  r\  direction 

Fv  viscous  flux  vector  in  the  rj  direction 

G  inviscid  flux  vector  in  £  direction 

Gv  viscous  flux  vector  in  the  ^  direction 

I  identity  matrix 

J  transformation  Jacobian 

k  thermal  conductivity 

L  model  length,  mm 

M  local  Mach  number 

p  pressure 

Pr  Prandtl  number,  0.72  for  air 

Q  flow  variable  vector 

q  heat  flux  vector 

R  gas  constant 

Re  Reynolds  number 

t  time 

T  absolute  temperature  (°K) 

u,v,w  Cartesian  velocity  components  in  the  x,y,and  z  directions. 

x,y,z  Cartesian  coordinates  in  streamwise,spanwise,  and  normal  directions 


a  angle  of  attack 

y  ratio  of  specific  heats 

p.  molecular  viscosity  coefficient 

q, rj>C  transformed  body  fitted  coordinates 

p  density 

x  shear  stress 

Subscripts 

aw  adiabatic  wall 

le  evaluated  at  the  wing  leading  edge 

n  normal  to  the  body  surface 

T  pitot  condition 

x  partial  derivative  with  respect  to  x 

y  partial  derivative  with  respect  to  y 

z  partial  derivative  with  respect  to  z 

%  partial  derivative  with  respect  to  £, 

t)  partial  derivative  with  respect  to  T) 

C,  partial  derivative  with  respect  to  £ 


AFIT/G  AE/AA/8  8  D-  34 


Abstract 


For  thin,  highly  swept  wings  operating  at  moderate  to  high  angles  of  attack,  the 
flow  over  the  wing  is  dominated  by  deformation  of  leading  edge  vortices.  These  vortices 
produce  a  minimum  pressure^hd  tftis  results  in  an  additional  lift  increment  •This  lift  in¬ 
crement  is  nonlinear  with  angle  of  attack  and  cannot  be  accurately  predicted  using  present 
design  methods.  ^ 

Y 

TheAhin-layer  Navier-Stokes  equations  were  used  to  calculate  the  flow  over  a 
straight  delta  wing  and  a  cranked  delta  wing.  The  straight  delta  wing  was  used  as  the  test 
case  due  to  the  availability  of  both  experimental  and  numerical  data.  Results  are  compared 
with  this  data  in  order  to  validate  the  numerical  procedure.  The  computer  code  uses  an  im¬ 
plicit,  time  marching  algorithm  developed  by  Beam  and  Warming.  The  solution  is  marched 
in  time  until  a  steady  state  is  achieved.  The  code  is  approximately  factored  and  diagonal¬ 
ized  in  order  to  reduce  computational  work.  A  solid  state  disk  is  used  in  order  to  allow  for 
the  large  grid  needed  for  a^dimensionaf  solution.  .. 

The  thin-layer  Navier-Stokes  equations  are  capable  of  accurately  calculating  vortical 
flows.  The  cranked  delta  wing  exhibited  flow  similar  to  a  straight  delta  wing  upstream  of 
the  crank.  The  vortex  generated  at  the  crank  quickly  became  paired  with  the  vortex  from 

the  front  of  the  wing.  Thevortex  location  aft  of  the  crank  changes  with  streamwise  loca- 

r.  . 

tion.  The  grid  resolution  is  important  when  trying  to  calculate  vortical  flows^  due  to  the 

Sy''  ■ 

large  gradients  in  both  tber  span  wise  and  normal  directions.  The  solid  state  disk  can  be 
used  to  run  problems  that  require  more  computer  memory  than  is  available.  Optimization 
of  the  program  input/output  should  be  done  for  running  the  code  with  the  solid  state  disk  in 

O'  ^ I 

drier  to  reduce  the  central  processor  unit  time  and  job  cost  1 

•  .  r\ 
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THIN-LAYER  NAVIER-STOKES  SOLUTIONS  FOR  A  CRANKED  DELTA  WING 


I  Imrodygmm 

The  current  class  of  high  performance  fighter  aircraft  must  be  able  to  operate  over  a 
wide  range  of  flight  conditions.  A  single  point  design  aircraft  is  no  longer  acceptable. 

They  must  have  efficient  supersonic  cruise  performance  and  yet  be  able  to  maneuver 
transonicly  for  dogfighting.  The  supersonic  requirement  has  driven  aircraft  designers  to 
use  thin,  highly  swept  wings.  Although  these  are  desirable  for  cruise  performance,  they 
greatly  degrade  the  maneuverability  of  the  aircraft.  This  performance  degradation  is  caused 
by  flow  separation  at  moderate  to  high  angles  of  attack.  This  difficulty  in  maintaining 
attached  flow,  has  led  designers  to  explore  the  possibility  of  using  vortex  lift  for  maneuver 
enhancement  (1:20). 

Current  fighters,  such  as  the  F- 16  and  the  F-18,  utilize  vortex  lift  to  enhance  their 
performance  by  the  use  of  wing-strake  planforms.  The  strakes  generate  strong  vortices, 
which  sweep  back  over  the  wing  and  generate  a  significant  increase  in  lift.  This  lift 
increment  is  a  result  of  the  decrease  in  pressure  on  the  upper  surface.  It  is  highly  nonlinear 
and  cannot  be  accurately  predicted  using  current  linear  design  methods. 

Newsome  and  Kandil  (2:2-3)  have  classified  the  flow  over  a  body  into  four  general 
categories.  The  first  category  is  for  attached  flow  which  occurs  at  low  angles  of  attack. 

The  second  category,  for  moderate  to  high  angles  of  attack,  is  identified  by  the  formation  of 
large  vortices  on  the  lee-side  of  delta  wings,  swept  wings,  low  aspect  ratio  wings,  and 
slender  bodies.  A  typical  flow  pattern  is  shown  in  Figure  1.  The  flow  is  dominated  by  the 
large  primary  vortices  generated  at  the  wing  leading  edge.  In  addition,  a  secondary  vortex 
is  formed  because  the  spanwise  flow  is  unable  to  negotiate  the  adverse  pressure  gradient. 
The  vortices  in  this  region  are  both  stable  and  symmetric  and  it  is  for  this  reason  that  most 
of  the  research  done  on  vortex  flows  has  been  in  this  region.  The  third  category  of  flow  is 


Figure  1 .  Vortex  Flow  Over  a  Delta  Wing  (3) 


for  very  high  angles  of  attack.  The  flow  in  this  region  produces  either  unstable  or 
asymmetric  vortices.  The  vortices  may  burst  or  become  asymmetric,  producing  large 
changes  in  the  body  forces  and  moments.  This  region  can  be  characterized  as  a  transition 
region  from  steady  stable  vortex  flows  to  unsteady,  asymmetric  flows.  The  majority  of  the 
research  in  this  area  is  related  to  predicting  vortex  breakdown  or  vortex  asymmetry.  The 
final  region  is  for  extreme  angles  of  attack.  The  flow  is  characterized  by  an  unsteady 
diffuse  wake  that  may  produce  vortex  shedding.  This  area  is  not  of  much  interest,  since 
most  flight  vehicles  are  not  likely  to  operate  in  it. 

For  sharp  delta  wings,  Stanbrook  and  Squire  (4)  have  shown  that  the  flow  can  be 
classified  according  to  the  normal  Mach  number  (Mn)  and  the  normal  angle  of  attack.(ctn), 
given  by 

Mn  =  Mo.V  l~sin2Acos2a  (1) 

an  =  tan-1(ianJ3Lj  (2) 

'cos  A' 

where  A  is  the  wing  leading  edge  sweep,  a  is  the  angle  of  attack,  and  Moo  is  the  freestream 
Mach  number.  The  Stanbrook-Squire  boundary  classifies  leading  edge  flows  into  either 
separated  flows  or  attached  flows  depending  upon  whether  the  normal  Mach  number  is  less 
than  or  greater  than  one. 

Miller  and  Wood  (5)  have  refined  this  boundary  and  identified  six  possible  types  of 
flow,  which  are  shown  in  Figure  2.  For  normal  Mach  numbers  less  than  one  and  small 
angles  of  attack,  the  flow  is  characterized  by  a  separation  bubble  at  the  leading  edge.  As 
angle  of  attack  is  increased  the  flow  separates  at  the  leading  edge  and  forms  a  vortex,  which 
is  often  referred  to  as  "classical  leading-edge  separation"  (2:3).  Normally  a  secondary 
vortex  will  form  under  the  the  primary  vortex.  This  is  due  to  the  spanwise  flow  being 
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Figure  2.  Possible  Types  of  Vortex  Flows  Over  Sharp  Delta  Wings  (5) 


unable  to  negotiate  the  adverse  pressure  gradient,  thus  separating  from  the  wing  surface. 

For  normal  Mach  numbers  greater  than  one,  there  are  three  types  of  flows  possible.  At  low 
angles  of  attack  the  flow  is  attached  and  the  leading  edge  expansion  is  terminated  by  a 
crossflow  shock.  With  increasing  angle  of  attack  the  shock  strengthens  and  causes  the 
formation  of  a  shock  induced  separation  bubble  in  board  on  the  delta  wing.  Further 
increase  in  angle  of  attack  causes  the  strengthing  of  a  thin  separation  bubble  with  a 
crossflow  shock  coalesced  on  top  of  it.  The  final  category  of  flow  occurs  for  very  high 
angles  of  attack,  with  either  subsonic  or  supersonic  normal  Mach  numbers.  This  flow  is 
characterized  by  strong  leading  edge  vortices  and  strong  crossflow  shocks. 

Methods  which  predict  vortical  flows  can  be  generally  classified  into  two 
categories.  The  first  category  models  the  vortex  in  an  approximate  manner  and  includes 
such  methods  as  the  Polhamus  suction  analogy,  vortex  lattice  methods  and  panel  methods. 
The  second  category  captures  the  vortex  as  a  solution  to  the  governing  equations  and 
includes  Euler  solutions  and  Navier-Stokes  solutions.  The  suction  analogy  (1,6)  is  an 
empirical  method  which  can  predict  forces  and  moments  but  is  unable  to  predict  pressures 
and  velocities.  Vortex  lattice  methods  (7-8)  are  able  to  predict  pressure  and  velocities  by 
modeling  the  vortex  sheet  as  discrete  line  vortices.  In  panel  methods  (9-13)  the  vortex 
sheet  is  more  accurately  modeled,  but  prior  knowledge  of  the  vortex  structure  is  required 
(14:825-826).  The  Euler  solutions  (15-19)  can  reproduce  the  primary  vortex  but  are  unable 
to  reproduce  the  secondary  structure.  Although  the  Navier-Stokes  solutions  (20-25) 
reproduce  all  facets  of  the  flow,  they  have  the  drawback  of  being  difficult  to  obtain  and 
computationally  expensive. 

Most  of  the  early  research  done  on  vortical  flows  used  one  of  the  methods  that 
approximately  modeled  the  vortex  core.  As  the  speed  and  size  of  computers  increased, 
Euler  and  Navier-Stokes  solutions  became  more  feasible.  The  Euler  equations  have  the 
advantage  of  being  simpler  to  solve  than  the  Navier-Stokes  equations.  The  major 
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disadvantage  to  using  the  Euler  equations  to  model  vortical  flows  is,  that  an  in  viscid 
approach  is  used  to  model  a  predominatly  viscous  phenomenon.  The  Euler  equations  can 
only  capture  the  primary  vortex  and  it  is  for  this  reason  that  the  Navier-Stokes  equations  are 
chosen  as  the  governing  equations  for  this  investigation. 

The  objective  of  this  research  is  to  use  the  Navier-Stokes  equations  to  calculate  the 
flow  over  a  cranked  delta  wing.  A  partial  listing  of  the  available  experimental  data  is 
contained  in  References  25-31.  The  configuration  chosen  corresponds  to  that  tested  by 
Henke  (31).  This  configuration  was  chosen  because  its  sharp  leading  edges  and  flat  upper 
surface,  will  produce  a  strong  vortex  which  is  free  of  body  influences.  The  test  conditions 
are  for  a  freestream  Mach  number  of  2.5  and  an  angle  of  attack  of  10°. 

The  solution  of  the  full  Navier-Stokes  equations  requires  a  substantial  amount  of 
computing  resources,  therefore  the  following  assumptions  were  made.  The  first  was  that 
the  viscous  effects  arc  only  significant  in  a  thin  layer  near  the  body.  The  viscous 
derivatives  normal  to  the  body  are  large  compared  to  the  viscous  derivatives  along  the 
body,  as  long  as  the  flow  is  attached  or  only  mildly  separated.  The  only  viscous 
derivatives  retained  are  those  normal  to  the  body  and  the  resulting  equations  are  known  as 
the  thin-layer  Navier-Stokes  equations  (32:9).  The  second  assumption  was  that  the  flow 
was  laminar.  Although  this  may  seem  restrictive,  the  laminar  numerical  results  of  Rizzetta 
and  Shang  (33)  and  Buter  and  Rizzetta  (34)  were  in  good  agreement  with  experimental 
data.  The  laminar  flow  assumption  appears  to  be  valid  for  freestream  Reynolds  numbers 
up  to  about  one  million. 

The  computer  code  used  is  the  ARC3D  code  written  by  Dr.  Thomas  Pulliam  of  the 
NASA  Ames  Research  Center  (32).  The  code  is  based  on  the  Beam  and  Warming  implicit 
approximate  factorization  algorithm  and  uses  the  thin-layer  approximation.  It  has  been 
modified  by  Dr.  Phil  Webster  of  the  Flight  Dynamics  Laboratory.  The  most  significant 
coding  change  has  been  the  incorporation  of  boundary  conditions  that  allow  for  a  branch 


cut,  between  the  upper  and  lower  body  surfaces.  This  allows  for  the  use  of  an  H-grid 
topology  which  gives  better  modeling  of  the  sharp  leading  edge. 

This  version  of  the  code  has  never  been  run  on  the  Cray  XMP- 1 2  computer  at 
Wright-Patterson  AFB.  Before  the  cranked  wing  calculations  could  be  performed,  the  code 
and  hardware  needed  to  be  checked.  This  was  done  by  running  a  test  case  for  which 
experimental  and  numerical  data  were  available.  The  test  case  corresponds  to  the 
configuration  tested  by  Monnerie  and  Werle  (35).  Available  numerical  data  includes  the 
full  Navier-Stokes  calculations  of  Rizzetta  and  Shang  (33)  and  Buter  and  Rizzetta(34),  and 
the  conical  solutions  of  Thomas  and  Newsome  (36)  and  Vigneron,  Rakich,  and  Tannehill 
(37). 


7 


H  Analysis 


Governing  Equations 

The  rime-dependent  compressible  Navier-Stokes  (NS)  equations  describe  the  con¬ 
servation  of  mass,  momentum,  and  energy  for  a  flowing  fluid.  Ignoring  body  forces  and 
external  heat  addition,  the  equations  can  be  written  in  non-dimensional,  strong  conserva¬ 
tion-law  form  as 
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The  equation  of  state  is  used  to  close  the  system  of  equations.  For  a  perfect  gas,  it  can  be 
expressed  in  terms  of  the  internal  energy  and  the  density  as 


p=(y-l)pei 


O' 


The  molecular  viscosity  (p.)  is  related  to  the  temperature  through  Sutherland's  formula  and 
has  the  units  of  kg/(m  s) 


p.  =  1.458(10) 


-6  T1S 
T-t-  110.4 


(11) 


The  preceding  equations  have  been  nondimensionalized  using  the  following  relations 
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where  the  dimensional  variables  are  denoted  by  an  asterisk,  freestream  conditions  by  <» 
and  L  is  the  body  reference  length.  Stokes's  hypothesis  is  used  to  relate  X.  to  p.  with  X. 
equal  to  -2/3  p.  (38:160-161). 

Coordinate  Transformation 

In  order  to  make  the  governing  equations  applicable  for  arbitrary  geometries,  they 
are  transformed  into  generalized  coordinates  (5,rj,C>t)  by  means  of  the  following  transfor¬ 
mation 


x  =  t 

5  =  5(x,y,z,t) 

T)  =  T|(x,y,z,t)  (13) 

C  =  C(x,y,z,t) 
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_LL_ 


Pulliam  (38:160-161)  has  shown  that  the  transformed  equations  written  in  strong  conser¬ 
vation-law  form  and  in  terms  of  the  contravaiiant  velocities  (U,V,W)  are  given  by 


5-jQ  +  +  9f|F  +  d^G  =  (9^EV  +  Bt|Fv  +  d^Gv) 


P 

^  Pu 
Q  =  J-!  pv 
pw 
e 


pU 

puU  +  ^xp 
E  =  J'1  pvU  +  £yp 
pwU  +  £zp 
_(e+p)U  +  ^tp. 


pv 

puV  +  T|xp 
F  =  J'  pvV  +  qyp 

pwV  +  TJzp 

_(e+p)V  +  Thp 


G  =  J1 


4x”^xx  ^jy^xy  +  ^sz^xz 
Fy  —  J  ^x^yx  "*■  ^y^yy  ^z^yz 
^x^zx  +  ^y^zy  +  &zz 

_  ^xPx  +  ^yPy  ■*"  ^zPz 


pW 

puw  +  Cxp 
PVW  +  CyP 
PwW  +  CzP 
_(e+p)W  +  ;iP_ 


fv  =  r1 


TJx'txx  +  ^ly’txy  +  ^Iz^xz 
Tlx^yx  +  'Hy'tyy  +  ^Iz^yz 
tlx^zx  +  Tly'Czy  +  TJz'tzz 
tlxPx  +  %Py  +  TlzPz 


0 


Cx'txx  "*■  Cy^xy  +  Cz^xz 
Cx^yx  +  Cy’fyy  Cz^yz 
Cx^zx  +  Cy^zy  +  Cz^zz 
CxPx  +  Cypy  +  CzPz  _ 


Px  =  7p.Pr_1  3xe,  +  utxx  +  vtxy  +  wtxz 
Py  =  TVlPr'1  9yei  +  UTyX  +  VTyy  +  wtyz 
Pz  =  Y)J.Pr_1  3zCj  +  u  ta  +  VTzy  +  wta 

u  =  +  ^XU  +  ^yV  +  CZW 

V  =Tlt  +  11xU+TlyV+TlzW 
W  =  Ct  +  CxU+CyV  +  CzW 


(18) 


(19) 


The  details  of  the  transformation  to  generalized  coordinates  is  contained  in  Appendix  A. 

Thin-Layer  Approximation 

The  thin-layer  approximation  applied  to  the  NS  equations  consists  of  neglecting  the 
viscous  derivatives  parallel  to  the  body  and  retaining  the  derivatives  normal  to  the  body. 
This  is  reasonable  for  high  Reynolds  number  flows,  as  long  as  the  flow  is  attached  or  only 
mildly  separated  (38: 160).  The  approximation  assumes  that  the  body  is  mapped  on  to  a 
Constant  plane,  and  the  viscous  derivatives  in  the  4  and  directions  are  neglected.  All  the 
convective  terms  are  retained  along  with  the  unsteady  terms.  Applying  the  thin-layer 
approximation  to  Equations  14-17,  Pulliam  and  Steger  have  shown  that  the  equations  can 
be  written  as 


9TQ  +■  e^E  +  9TjF  +  9^G  — 


1 

—  3rS 
Re  ? 


(20) 
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where  Q,  E,  F,  and  G  are  the  same  as  before.  The  viscous  flux  terms  parallel  to  the  body 
are  neglected  (3^EV  =  d-^Fy  =  0)  and  the  viscous  flux  term  normal  to  the  body  is  given  by 


0 


S  =  j-1 


lanuu^  +  (M-/3)m2Cx 
^miv?  +  (d/3)m2Cy 
M-miw^  +  4t/3)m2£2 
}imim3  +  (|x/3)  m2nu 


m,=d+£+;« 

m2  =  +  Cy  +  Czw; 

m3  =  ~3^(u2  +  v2  +  w2)  +  Y 

m4  =  CxU  +  CyV  +  Czw 


(21) 


(22) 


The  thin-layer  Navier-Stokes  (TLNS)  equations  are  a  mixed  set  of  hyperbolic- 
parabolic  partial  differential  equations  (PDE)  in  time.  These  have  the  same  form  as  the  NS 
equations,  and  can  be  solved  using  the  same  type  of  methods.  If  the  unsteady  terms  are 
dropped,  the  equations  become  a  mixed  set  of  hyperbolic-elliptic  equations.  These  equa¬ 
tions  are  more  difficult  to  solve  and  therefore  most  solutions  for  the  compressible  NS 
equations  have  used  the  unsteady  form.  The  steady-state  solution  is  obtained  by  marching 
the  solution  in  time  until  converged,  and  is  known  as  the  time-dependent  approach.  The 
time-dependent  approach  can  be  explicit  or  implicit  and  is  usually  second-order  accurate  in 
space.  If  time  accuracy  is  desired,  the  algorithm  should  be  second-order  accurate  in  time. 
For  steady-state  calculations,  a  first-order  accurate  algorithm  can  be  used  to  accelerate  con¬ 
vergence  (39:424,482). 
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The  pressure  at  the  body  surface  can  be  found  using  a  normal  pressure  momentum 
relation.  Pulliam  has  shown  (38:161)  by  combining  the  three  transformed  momentum 
equations,  the  normal  pressure  is  given  by 

p„(d  +  d  +  d)  -  p(a&  +  +  va^y  +  wa^J  - 

pir(C„U5  +  CyV,  +  tywj  -  pV(;xU^  +  tyv„  +  Cywj  (23) 

and  the  Cartesian  velocities  can  be  determined  from 

'  u  '  (%£z  ~  HzCy)  — (^yCz  ~  ^zCy)  (^Hz  ~  ^z^y)  U  - 

v  -I’1  -kC.-1.tJ  («.-«»)  -(titl.-t.1lJ  V_n,  (24) 

LwJ  [  (n.Cy-tiyC.)  -(t»C,-WJ  (t.iy-t,n.)  J[  . 

The  above  expressions  are  valid  for  inviscid  or  viscous  flows  and  for  steady  or  unsteady 
body  motion. 

The  no  slip  boundary  condition  is  specified  at  the  body  surface,  where  U=V=W=0 
and  u,v,  and  w  are  determined  from  Equation  24.  In  the  case  of  a  stationary  grid  where 
c,t=rii=^t=0,  u,v,  and  w  are  zero.  The  surface  pressure  is  calculated  by  integrating 
Equation  23.  For  viscous  flow  over  a  stationary  grid,  Equation  23  reduces  to  the  pressure 
gradient  normal  to  the  surface  being  zero  (Pn=0).  The  body  is  assumed  to  be  adiabatic 
(Tn=0)  and  the  density  is  found  by  linear  extrapolation  from  the  point  above.  The  farfield 
boundary  conditions  are  specified  to  be  frees tream  values,  except  for  the  downstream 
condition.  A  first-order  extrapolation  was  used  as  the  downstream  boundary  condition.  A 
plane  of  symmetry  was  imposed  at  the  centerline.of  the  body,  where 


i* - Cf.-Tk  . 


Ill  Numerical  Procedure 

The  numerical  procedure  used  to  obtain  solutions  of  the  governing  equations  is  out¬ 
lined.  Topics  include  the  choice  of  algorithm,  and  formulation  of  the  finite-difference  equa¬ 
tions.  Emphasis  is  placed  on  detailing  the  simplifications  used  (approximate  factorization 
and  diagonalization)  and  their  impact  on  code  efficiency  and  accuracy.  The  implementation 
of  the  nonlinear  dissipation  model  is  discussed.  The  solution  process  for  one  time  step  is 
outlined  for  running  the  code,  in  core  or  out  of  core  with  the  Solid-State  Disk  (SSD). 


Implicit  Time  Matching  Algorithm 

A  time  marching  finite-difference  scheme,  developed  by  Beam  and  Warming 
(40:1 18-120),  is  used  to  solve  the  thin-layer  Navier-Stokes  equations.  This  is  an  alter¬ 
nating-direction  implicit  (ADI)  scheme  and  is  similar  to  schemes  developed  by  Lindemuth 
and  Killeen  (41)  and  McDonald  and  Briley  (42).  The  scheme  uses  an  implicit,  three-point, 
time-differencing  formula  in  the  form 


(26) 


the  parameters  and  cp  can  be  chosen  to  produce  a  scheme  which  is  either  first-  or  second- 
order  accurate  in  time.  Since  we  are  primarily  interested  in  steady-state  solutions,  a  first- 
order  scheme  is  chosen.  If  d=l  and  cp=0,  this  results  in  the  Euler  implicit  scheme  (39:490) 

AQn  =  At^(AQn)  +  At^(Qn)  +o[(At)2]  (27) 


^  n  n  +  l  ^  || 

and  AQ  =  Q  -  Q  ,  thus  Equation  27  can  be  rewritten  as 


AQ"-At|(Q"*')=0 


(28) 


lb 


Writing  Equation  20  at  the  n+1  time  level  yields 


34"l=-a5E"'-311F"*1-a5G"*1 


-L3rS 

Re  ^ 


n  +  1 


(29) 


Finally,  substitution  of  Equation  29  into  Equation  28  yields 

AQn  +AtJa^En+1  +a^Fn+1  +8;Gn+1-^-a^Sn+1)  =  0  (30) 


Linearization 

For  Equation  30  to  be  solved  for  Q  the  flux  vectors,  E,  F,  G,  and  S,  which  are 
nonlinear  functions  of  Q,  must  be  linearized.  The  inviscid  flux  vectors  can  be  linearized  in 
time  by  using  a  Taylor  series  about  Qn  (32:12): 

E"*1  =  g"  +  ^  AQn  =  £”  +  A"  AQ" 

3Qn 

Fn+1  =  Fn  +  ^-AQ"  =  Fn  +  BnAQn  (31) 

aQn 

Gn+1  =  Gn  +  AQn  =  Gn  +  CnAQn 

aQn 

Steger  has  shown  (43:3-4)  that  the  viscous  flux  vector  can  be  linearized  using  a  Taylor 
series  and  is  given  by 

Sn+l  =  Sn  +  J- 1  ^|— J  AQn  =  Sn  +  J-1  MnJAQn  (32) 

aQn 
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The  flux  Jacobians  arc  given  in  Appendix  B.  Substituting  Equations  31  and  32  into 
Equation  30  yields  the  "delta  form"  of  the  algorithm 


[i  +At3^An  +  At9r,Bn  +  Ata^C"  -  Re"1  At9;J-1Mnj]  AQn  = 

-At(a^En  +  a^F”  +  a^Gn  -  Re-1  acS  ")  (33) 

In  Equation  33  that  the  left-hand-side  contains  the  unknown  AQ  and  is  sometimes 
referred  to  as  the  "implicit"  part.  The  right-hand-side  contains  the  known  quantities  and  is 
referred  to  as  the  "explicit"  pan  (32:13).  Applying  second-order,  central-difference  oper¬ 
ators  to  equation  33  yields  the  final  form  of  the  time  marching  algorithm: 

[l  +At5^An  +  AtS^B"  +  AtS^C"  -  Re"1  AtS^-MM".!]  AQn  = 

-At(8^En  +  5T1Fn  +  8;Gn-Re-15cSn)  (34) 


Approximate  factorization 

To  integrate  the  full  three-dimensional  operator  would  be  prohibitively  expensive. 

1  One  simplification  that  can  be  used  is  to  approximately  factor  the  three-dimensional  oper¬ 

ator  into  three  one-dimensional  operators.  If  terms  on  the  0(At2)  are  neglected  the  left- 

hand-side  of  Equation  34  is  factored  into  (32:79) 
i 

[i  +  AtS^A"]  [i  +  AtS^B"]  [i  +  AtS^C"  -  Re'1At8<;J-1Mnj]AQn  (35) 

1  Neglecting  the  0(At2)  terms  will  not  degrade  the  accuracy  of  the  scheme  since  it  is  first 

order  accurate  in  time.  Equation  34  can  now  be  written  using  approximate  factorization  as 


[l-)-At8^An][l  +  At8?1Bn][l  +  At8?Cn-Re-1At8?J-1Mnj]AQn  = 


-At(5^En  +  S^F"  +  8^G"  -  Re"1  8^Sn) 


(36) 


The  approximate  factorization  reduces  the  (Jmax-KmaxLmax-5)  x  (Jmax-Kmax  Lmax-5) 
banded  matrix  down  to  a  set  of  three  block  tridiagonal  matrices.  The  size  of  any  matrix  is 
now  at  most  (max  I  Jmax,Kmax,Lmax  1-5)  x  (max  I  Jmax,Kmax,Lmax  1-5).  The  solution 
now  consists  of  three  sweeps,  one  in  the  £  direction,  one  in  the  T(  direction,  and  one  in  the 
£  direction.  Each  block  tridiagonal  matrix  can  be  solved  using  a  block  lower-upper  decom¬ 
position  (LUD)  (32:17). 

Pia&pnalizflpQn 

Approximate  factorization  was  used  to  reduce  the  block  three  dimensional  implicit 
operator  of  Equation  34  down  to  three  one-dimensional  block  tridiagonal  matrices.  The 
solution  of  these  matrices  are  still  computationally  complex,  since  they  involve  the  solution 
of  5x5  blocks.  One  way  to  decrease  the  work  is  to  decouple  the  equations.  If  the  operators 
are  diagonalized,  the  block  structure  is  decomposed  into  five  scalar  operators. 

Diagonalization  of  the  Euler  equations  is  presented  and  then  the  method  is  extended 
to  the  Navier-Stokes  equations.  This  is  due  to  the  fact  that "  the  viscous  flux  Jacobian  M" 
is  not  simultaneously  diagonalizable  with  the  flux  Jacobian  Cn  (32:19)."  Warming,  Beam, 
and  Hyett  (44:1037-1041)  have  shown  that  the  inviscid  flux  Jacobians  can  be  diagonalized 
since  they  have  real  eigenvalues  and  a  complete  set  of  eigenvectors.  The  flux  Jacobians  are 
written  in  terms  of  their  eigenvalues  and  eigenvectors  as 

Bn*(TT,A„T^)n  (37) 
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te&fe  <KU  »> 


where  is  the  eigenvector  matrix  of  A  and  likewise,  Tn  for  B,  and  Tj;  for  C.  The  eigen¬ 
vector  and  eigenvalue  matrices  are  written  out  in  Appendix  C.  Equation  36  can  now  be 
written  neglecting  the  viscous  flux  Jacobian  as 

+a^t5a5t-1F][(t„t;1F  +  a^t,a,t;'F] 

x  [(t^T'1)"  +  A  t5<|T  7 1 1*]  AQn  =  Rn  (38) 

where 

Rn  =  -At(8^En  r  S^F"  +  8^Gn  -  Re-^S")  (39) 

The  eigenvector  matrices  are  factored  outside  of  the  difference  operators  and  this  yields  the 
"diagonal"  form  of  the  algorithm. 

Tjl  +  At6^A^]Nn[l  +  At5nAj|]?n[l  +  At6<;A  f  AQn  =  Rn  (40) 


where 


Pn=(T^1Tr)n  (41) 

Since  the  eigenvector  matrices  are  functions  of  £,r|,  and  factoring  them  outside 
of  the  operator  introduces  an  error.  Pulliam  and  Chaussee  (45:356-359)  have  shown  that 
for  steady-state  solutions,  the  accuracy  of  the  solution  is  not  affected.  The  reason  being, 
for  steady-state  solutions  AQ  approaches  zero.  Since  the  right-hand-side  is  unchanged  by 
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the  diagonalization  the  error  is  determined  by  the  order  of  the  differencing  chosen.  They 
have  also  shown  the  time  accuracy  of  the  algorithm  to  be  at  most  first-order. 

The  above  discussion  only  applies  to  the  Euler  equations  since  the  viscous  flux 
Jacobian  has  been  neglected.  To  extend  diagonalization  to  the  Navier- Stokes  equations, 
Pulliam  suggests  (32:20)  including  a  diagonal  term  on  the  implicit  side  to  approximate  the 
viscous  Jacobian  eigenvalues.  The  current  estimates  are 

Xvk)  =  p|iRe-1J-I(4x  +  +  Sz) 

Xv(Tl)  =  ppRe- 1 J'  +  Til)  (42) 

^v(c)  =  PHRe-1J-1(Cx  +  Cy  +  Cz) 

For  the  thin-layer  approximation  the  ^  and  T|  terms  arc  ignored  and  the  Xv(Q  term  is  added 
to  the  eigenvalues. 

AaifidalPissipaQoii 

By  the  use  of  linear  stability  analysis,  it  can  be  shown  that  the  diagonalized  algo¬ 
rithm,  Equation  40,  is  unconditionally  stable.  In  reality  this  is  not  true,  especially  when  the 
system  is  nonlinear.  "Scales  of  motion  appear  which  cannot  be  resolved  by  the  numerics 
and  are  due  to  the  nonlinear  interactions  in  the  convection  terms  of  the  momentum 
equations  (32:26)." 

One  way  of  dealing  with  these  numerical  instabilities  is  to  add  some  numerical 
dissipation  to  the  algorithm  which  does  not  alter  the  accuracy  of  the  solution.  The  dis¬ 
sipation  model  chosen  is  a  mixed  second  order,  fourth  order  model  due  to  Jameson  (46). 
The  nonlinear  model  for  the  £  coordinate  is  given  by 

Dr  =  ♦  aiuJ-uK A  -  <$AVA)J  <«> 
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where 


£j 5.1  *  «C2Atm«x(Ti.k4+i,T j.k.i»Y j.k.1-1) 

r.  -  ipj.U^i  ~  2Pj,k,i  +  Pj.k.1-1  j 
j,k’1  |Pj,k,l+l  +  2Pj.kJ  +  Pj.U-l| 

e(j4k.i  =  max(°-^At-ej2u) 

0j,k,l  =  I  W|  +  a(^x  +  Cy  +  Cz) 


(44) 


and  a  is  the  spectral  radius  scaling  for  C.  and  are  the  first  order  accurate  backward 
and  forward  difference  operators.  Typical  values  of  the  constants  K2  and  K4  are  .25  and 
.01,  respectively.  Similar  terms  are  used  for  the  £  and  T|  directions.  Applying  this  model 
to  both  the  implicit  and  explicit  sides  of  Equadon  40  yields 

Tj[l  +  AtS^A?  -  AtDjjjN^I  +  AtSrjA^  -  AtDnl]p"[l  +  AtS^A;  -  AtD;l](f^1)n AQn  = 
-AtjS^E"  +  S^F"  +  5?Gn  -  Re'1 5?S")  -  (D^  +  D^  +  D^)Qn  (45) 

Since  the  dissipation  is  a  fourth  order  model,  this  necessitates  the  use  of  scalar  penta- 
diagonal  solvers.  Although  pentadiagonal  equations  are  more  complicated  to  solve  than 
tridiagonal  equations,  the  use  of  implicit  dissipation  improves  convergence  and  "robust¬ 
ness"  of  the  algorithm  (47:16). 

Solution  Procedure 

Since  the  algorithm  makes  use  of  approximate  factorization,  the  solution  can  be 
obtained  by  sweeping  in  the  ^,ri,  and  £  directions  sequentially.  The  solution  process 
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becomes  a  series  of  matrix-vector  multiplies  and  solutions  of  scalar  pentadiagonal  equa- 
tions.This  process  consists  of  the  following  eight  steps  at  each  time  level. 

1 .  A  matrix-vector  multiplication  at  each  grid  point 

MT5T[ft”-(<VDn  +  D?£”] 

2.  The  solution  of  five  scalar  pentadiagonal  equations  for  the  ^  direction. 

§2  =  [l  +  AtS^Ajs  -  AtD^l]-1  S\ 

3 .  A  second  matrix-vector  multiplication  at  each  grid  point 

S3=(n  ‘I  s2 

4.  The  solution  of  five  pentadiagonal  equations, this  time  in  the  r\  direction. 

S4  =  [i  +  AtSqA^  -  AtD^l]"1  S3 

5 .  A  third  matrix-vector  multiplication  at  each  grid  point 

a  /a,  i  \n 

Ss  =  (P  ‘I  S4 

6.  The  solution  of  five  scalar  pentadiagonal  equations,  in  the  £  direction. 

S6  =  [l  +  At5^-AtD;l]_1S5 

7 .  The  final  matrix -vector  multiplication  at  each  grid  point. 

AQn  =  (T?)"S6 

8 .  The  solution  is  then  updated  by 

Q"*1  =  Qn  +AQ" 
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Solid  State  Disk  Implementation 

The  calculations  were  done  on  a  Cray  XMP-12  computer  using  a  Solid  State  Disk 
(SSD).  The  SSD  was  used  because  of  the  large  grid  and  large  amounts  of  storage  required 
for  the  solution  of  a  three  dimensional  problem. 

The  solution  process  using  the  SSD,  entails  two  sweeps  of  the  domain  for  each 
time  step.  The  first  sweep  is  for  planes  of  constant  r\  (k  planes),  where  the  differencing  in 
the  C,  direction  is  done.  The  second  sweep  is  for  planes  of  constant  £  Q  planes),  where  the 
4  and  rj  differencing  is  done. 

The  SSD  is  set  up  with  two  working  files.  The  first  file  contains  the  metrics  and 
the  second  file  contains  the  flow  variables.  At  each  plane  the  metrics  and  the  flow  variables 
are  read  into  the  in-core  memory.  The  calculations  are  carried  out  and  the  updated  flow 
variables  are  loaded  back  on  to  the  SSD.  There  is  no  need  to  unload  the  metrics  since  they 
are  invarient  with  time.  The  algorithm  advances  to  the  next  plane  and  the  process  is 
repeated. 

Using  the  SSD  is  one  way  of  running  problems  which  require  more  core  memory 
than  may  be  available.  The  trade  off  is  that  the  code  is  not  as  efficient  as  a  code  that  runs 
totally  in  core.  Appendix  D  contains  a  comparison  of  Central  Processor  Unit  (CPU)  times, 
input/output  (I/O),  and  storage  required  to  run  the  code  in  core  and  out  of  core  using  the 
SSD. 


IV  Results  and  Discussion 

The  diagonalized  algorithm  is  used  to  calculate  the  flow  over  a  75°  delta  wing  and  a 
80°/69°  cranked  delta  wing  at  angle  of  attack.  The  delta  wing  calculations  are  compared 
with  the  experimental  data  of  Monnerie  and  Werle  (35)  and  the  numerical  results  of  Buter 
and  Rizzetta  (34).  The  cranked  wing  solution  is  compared  with  experimental  data  of  Henke 
(31).  The  calculations  were  done  on  the  Cray  XMP- 12  computer  at  Wright-Patterson  Air 
Force  Base  (AFB),  Ohio.  The  information  and  files  required  to  run  this  code  at  the  Air 
Force  Institute  of  Technology,  are  contained  in  Appendix  E. 

Delta  Wine 

Since  this  version  of  the  code  had  not  been  previously  run  using  the  Cray  at 
Wright-Patterson  AFB,  it  was  necessary  to  insure  there  w ere  no  hardware  oar  system 
software  incompatibilities.  A  75°  delta  wing  at  a  fneestream  Mach  number  of  1.95  and 
angle  of  attack  of  10°  was  chosen  as  the  test  case.  This  configuration  was  tested  by 
Monnerie  and  Werie  (35)  and  the  model  geometry  is  given  in  Figure  3.  Numerical 
solutions  using  the  full  Navier-Stokes  equations  and  the  conical  approximation  to  the 
Navier-Stokes  equations  are  contained  in  references  31-34.  Comparisons  will  be  made 
with  the  experimental  data  and  with  the  numerical  results  of  Buter  and  Rizzetta  (34). 

Due  to  the  simple  geometry  of  the  delta  wing,  the  grid  was  generated  using  an 
algebraic  grid  generator.  The  grid  consisted  of  30  points  in  the  stream  wise  (%)  direction, 

55  points  in  the  span  wise  (tj)  direction  and  65  points  in  the  normal  (Q  direction.  A  typical 
stream  wise  grid  plane  (q,Q  is  shown  in  Figure  4.  The  grid  incorporates  a  branch  cut 
between  the  upper  and  lower  body  surfaces  to  allow  for  the  use  of  an  H-grid.  This 
approach  was  used  to  eliminate  the  small  radius  of  curvature  of  the  grid,  at  the  wingtip,  that 
results  from  using  a  C  grid.  Since  the  freestream  Mach  number  is  supersonic,  the  upstream 
propagation  of  disturbances  is  extremely  limited.  This  implies  that  there  is  no  requirement 


for  a  wake  region  and  the  domain  ends  at  the  wing  trailing  edge.  Orthogonality  is  not 
enforced  at  the  boundaries,  but  previous  work  (33,34)  has  shown  that  good  results  can  be 
obtained  using  this  approach.  Points  were  clustered  in  the  high  gradient  areas  such  as,  the 
nose,  the  wingtip  and  at  the  body  surface,  by  the  use  of  exponential  functions. 

Figure  5  shows  the  comparison  of  Pitot  pressures  with  the  experimental  data.The 
Pitot  pressures  have  been  normalized  by  the  freestream  Pitot  pressure.  The  thin-layer 
Navier-Stokes  calculations  show  good  agreement  on  the  location  and  the  strength  of  the 
primary  vortex.  The  comparison  with  the  work  of  Buter  and  Rizzetta  are  shown  in  Figure 
6.  Again  good  agreement  is  seen  with  the  full  Navier-Stokes  calculations.  As  was  seen  in 
the  previous  comparison,  the  thin-layer  calculations  exhibit  more  total  pressure  loss,  but  the 
crossplane  velocities  and  the  wing  surface  pressures  agree  well.  The  present  study  exhibits 
a  little  less  development  in  the  secondary  separation  region  than  do  the  full  Navier-Stokes 
solutions  as  seen  in  the  surface  pressure  plot 

The  streamwise  development  of  the  leading  edge  vortex  is  given  in  Figure  7. 
Traveling  from  the  nose  to  the  trailing  edge,  the  vortex  is  seen  to  strengthen.  This  is 
marked  by  a  decrease  in  pressure  of  the  primary  vortex  core  and  an  increase  in  the 
crossplane  velocities.  Also  from  the  crossplane  velocities,  the  formation  and  strengthening 
of  the  secondary  separation  is  seen,  although  this  is  not  significant  until  aft  of  60  percent  of 
the  wing  length.  Upon  closer  examination  of  the  crossplane  velocities,  it  is  noted  that  the 
center  of  the  vortex  and  gross  features  of  the  flow  do  not  change  rapidly  in  the  streamwise 
direction.  Although  the  flow  is  not  truly  conical,  it  is  changing  slowly  enough  that 
ignoring  the  viscous  effect  will  not  significantly  affect  the  overall  flow  field  structure.  The 
same  trend  can  be  seen  in  the  surface  pressure  plots  of  Figure  8.  The  primary  vortex  is 
seen  to  increase  in  strength  with  aft  movement. 
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Figure  6.  (corn)  Crossplane  Velocities  at  x/L=0. 8 


Full  Navier-Stokes 


Figure  6.  (corn)  Surface  Pressures  at  x/L=0.8 
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Figure  8.  Delta  Wing  Surface  Pressures 

x/l=0.2  (upper),  x/l=0.4  (lower) 


The  configuration  chosen  for  this  study  was  that  tested  by  Henke  (31)  and  is  shown 
in  Figure  9.  It  has  an  80°  leading  edge  sweep  forward  of  the  crank  and  is  followed  by  a 
69°  sweep.  The  crank  occurs  at  40  percent  of  the  body  length.  The  wing  has  a  triangular 
cross  section  with  a  flat  upper  surface  and  sharp  leading  edges.  The  angle  between  the 
upper  and  lower  surfaces  is  40°  and  this  results  in  a  relatively  thick  body  at  the  trailing 
edge. 

An  algebraic  grid  generator  was  also  used  to  fit  a  grid  around  the  cranked  delta 
wing.  The  grid  generator  was  originally  written  to  generate  grids  for  delta  wings  and  has 
been  modified  by  the  author  to  generate  grids  for  cranked  delta  wings.  Appendix  F 
contains  the  listing  of  the  grid  generator  used  for  these  calculations.  The  grid  consists  of 
127,400  points,  28  in  the  streamwise  direction,  65  in  the  span  wise  direction,  and  70 
normal  to  the  body.  A  typical  streamwise  grid  plane  is  shown  in  Figure  10.  As  can  be 
seen  from  the  figure,  points  have  been  clustered  at  the  body  surface  and  the  wing  tip.  The 
grid  consists  of  41  points  on  the  wing  surface  in  the  spanwise  direction  and  25  points  in  the 
streamwise  direction.  In  addition  to  clustering  at  the  nose,  points  were  clustered  in  the 
vicinity  of  the  wing  crank.  This  was  necessary  because  of  the  large  changes  in  the  flow, 
resulting  from  the  start  up  of  the  crank  vortex.  As  before  a  branch  cut  is  used  to  allow  for 
an  H  grid  topology.  Since  the  freestream  Mach  number  is  supersonic,  there  is  no  signif¬ 
icant  signal  propagation  from  the  wake  region,  therefore  the  domain  ends  at  the  wing 
trailing  edge. 

The  first  plane  of  the  grid  on  the  body  is  located  at  x/l=0.05  from  the  apex  of  the 
delta  wing.  It  was  found  that  placing  a  plane  of  data  farther  forward  than  this  caused  the 
program  to  crash.  It  is  thought  that  the  points  become  compressed  in  the  t|-£  plane  and  this 
causes  large  unsmooth  variations  in  the  metrics.  The  correction  to  this  problem  was  to 
move  the  initial  plane  rearward  and  open  up  the  r|-£  plane. 
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The  configuration  chosen  was  for  a  ffeestream  Mach  number  of  2.5  and  an  angle  of 
attack  of  10  degrees.  This  angle  of  attack  was  chosen  because  it  was  thought  to  be  high 
enough  to  develop  stable  vortical  flow.  The  ffeestream  Reynolds  number  is  686,000  and 
the  wind  tunnel  stagnation  conditions  are  10  psia  and  20°  C. 

Figure  1 1  shows  the  development  of  the  Pitot  pressures  and  cross  plane  velocities 
along  the  wing.  The  flow  over  the  wing  in  front  of  the  crank  is  similar  to  that  found  on 
straight  delta  wings.  From  the  crossplane  velocity  plot,  for  a  streamwise  location  of 
x/l=0.398,  it  is  seen  that  the  center  of  the  primary  vortex  is  located  at  approximately  50 
percent  of  the  wing  span.  The  location  of  the  vortex  center  remains  fairly  constant  with 
streamwise  location.  This  was  also  true  for  the  flow  over  the  75°  delta  wing  that  was 
shown  previously.  In  addition  to  the  primary  vortex  structure,  there  was  also  a  secondary 
structure  located  at  about  85  percent  of  the  wing  span. 

Moving  aft  of  the  crank,  a  second  vortex  forms  at  the  wing  crank.  This  second 
vortex  can  be  seen  in  the  plots  of  Figure  1 1,  at  a  wing  location  of  x/L =0.414.  The  vortex 
center  is  located  at  approximately  85  percent  of  the  wing  span.  The  area  of  secondary  sep¬ 
aration  seen  in  front  of  the  crank,  has  started  to  diminish.  Moving  aft  on  the  wing 
(x/L=0.453)  shows  no  evidence  of  the  second  vortex.  The  vortex  generated  at  the  crank 
appears  to  become  paired  with  the  first  vortex  and  the  secondary  vortex  has  now  com¬ 
pletely  disappeared  as  a  seperate  identity.  Finally,  looking  at  a  plane  well  aft  of  the  crank  , 
x/L=0.88,  an  area  of  low  Pitot  pressure  has  formed  between  the  first  vortex  and  the  wing 
leading  edge. 

The  vortex  location  is  seen  not  to  change  with  streamwise  location  in  front  of  the 
crank  but  this  is  not  true  aft  of  the  crank.  Figure  12  shows  a  comparison  between  the  flow 
in  front  of  the  crank  and  well  downstream  of  the  crank.  Aft  of  the  crank  the  vortex  center 
is  seen  to  move  inboard  and  towards  the  wing  surface.  The  area  of  secondary  separation 
present  in  the  flow  in  front  of  the  crank  has  been  replaced  by  an  area  of  low  Pitot  pressure. 
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The  surface  pressures  show  that  the  overall  magnitude  of  the  pressure  is  approximately  the 
same  for  both  planes.  The  primary  difference  is  that  the  flow  over  the  cranked  portion 
must  negotiate  a  larger  pressure  gradient  than  the  flow  in  front  of  the  crank. 

Figure  13  shows  a  comparison  of  the  upper  surface  pressures  with  experimental 
data.  The  upper  surface  pressure  exhibits  the  same  trend  as  the  experimental  data  but  is 
about  10  percent  less  for  the  outboard  wing  section.  It  appears  that  the  grid  may  be  too 
coarse  in  the  spanwise  and  normal  directions  to  adequately  resolve  all  the  details  of  the  flow 


Figure  1 1.  Pitot  Pressures  and  Crossplane  Velocities  for 
the  Cranked  Delta  Wing,  x/L=0.398 
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Figure  12.  (Cont)  Surface 


Experiment 


Cranked  Wing  Surface  Pressures  Compared  with  Experiment 


The  Thin-layer  Navier-Stokes  equations  were  shown  to  be  capable  of  reproducing 
both  the  primary  and  secondary  structures  of  vortical  flows.  For  flow  over  a  straight  delta 
wing,  good  agreement  was  obtained  with  the  full  Navier-Stokes  solutions  of  Buter  and 
Rizzetta.  The  present  study  did  exhibit  slightly  more  Pitot  pressure  loss  than  either  the 
experimental  data  or  the  full  Navier-Stokes  solutions. 

The  cranked  delta  wing  exhibited  flow  similar  to  a  straight  delta  wing  for  locations 
upstream  of  the  crank.  Vortex  location  did  not  vary  greatly  with  streamwise  location.  A 
second  vortex  formed  at  the  crank  and  then  quickly  became  paired  with  the  first  vortex. 
Moving  downstream  the  vortex  center  moves  inboard  and  towards  the  wing  surface.  The 
secondary  separation  present  upstream  of  the  crank  has  been  replaced  by  an  area  of  low 
Pitot  pressure.  The  magnitude  of  the  wing  surface  pressures  are  approximately  the  same 
for  the  flow  fore  and  aft  of  the  crank.  The  major  difference  is  the  flow  downstream  of  the 
crank  must  negotiate  a  larger  pressure  gradient  than  the  flow  upstream  of  the  crank. 

Grid  resolution  appears  to  be  critical  when  trying  to  calculate  vortical  flows.  It  is 
necessary  to  cluster  points  not  only  normal  to  the  body  but  also  in  the  span  wise  direction 
near  the  leading  edge.  This  is  required  in  order  to  accurately  resolve  the  large  gradients  in 
the  T)  direction.  The  result  of  insufficient  numerical  resolution  is  degraded  definition  of  the 
secondary  vortical  structure  or  even  a  failure  to  resolve  the  primary  structure.  For  a  delta 
wing  the  grid  can  be  fairly  coarse  in  the  streamwise  direction  as  long  as  the  vortices  are 
stable.  If  the  flow  becomes  unsteady,  such  as  with  vortex  bursting,  then  the  resolution  in 
the  streamwise  direction  should  be  fairly  fine.  The  above  is  only  true  when  there  is  little  or 
no  trailing  edge  effects  propagating  upstream.  The  cranked  delta  wing  will  require  a  fine 
grid  in  the  streamwise  direction  due  to  the  flow  changing  downstream  of  the  crank. 

The  use  of  the  Solid  State  Disk  (SSD)  allows  for  the  use  of  grids  that  require  more 
computer  memory  than  is  available.  The  accuracy  of  the  solution  is  not  affected  when 


using  the  SSD.  The  main  drawback  to  using  the  SSD  with  the  ARC3D  code  seems  to  be 
the  amount  of  input/output  required  If  the  input/output  can  be  better  optimized  this  would 
reduce  the  required  central  processor  unit  time  and  also  reduce  the  job  cost. 
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The  governing  equations  can  be  transformed  from  Cartesian  coordinates  (x,y,z,t) 
into  generalized  coordinates  (4/n,C,t)  by  the  following  transformation 


^  =  4(x,y,z,t) 


ti  =n(x,y,z,t) 
C  =  C(x,y,2,t) 


Using  the  chain  rule,  the  Cartesian  partial  differentials  can  be  written  in  terms  of  the 
transformed  coordinates 


d  _dx  d  3£  8  3n  3  3£  8  ^ 

- ■* - - *  =  3T  +  +  ritSp  +  Ldc 

3t  3t  3x  3t  32;  3t  dq  5t  8£ 

8  8x8  8^8  8ri8  8C8r. 

T-  =  7'~  +  ~~+— ~+7-— =  4x3^  +  Tlx^  + 

8x  8x  8x  8x  32;  8x  8q  8x  3£ 

8  8t8  32;  8  8r]8  3£  3  c- 

—  =  - - +  —  —  h — - —  +  ——  =  <*  3? +  1^3.  +  rv5c 

8y  8y  8x  8y  32;  8y  dq  3y  8£ 

3_8x3  82;  8  3t}3  3£  3 

3z  3z  8t  +  3z  32;  +  3z  3q  +  3C  ”  ^  ^  +  ^  ^  +  ^  C 


Applying  Equation  A2  to  Equation  3  yields 


3tQ  +  3^tQ  +  2;XE  +  ^yF  +  2;zg)  +  +  t|xE  +  qyF  +  t|zG)  + 

^(CtQ  +  CxE  +  CyF  +  CzG)  =  Re-I[8^xEv  +  2;yFv  +  ^ZGV)  + 

^TlxEv  +  TiyFv  +  T)ZGV)  +  3^XEV  +  £yFv  +  CzGv)  (A3) 
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the  metrics  appearing  in  the  above  equation  are  defined  as 

^  =  j(y^-y^) 

^=J(Z11X;-XT1Z;) 

^  =  j(xT,y;-yr1xJ 
Tlx  =  J  (z§y?  -  y&) 
qy  =  j(x^-x^) 

Tlz  =  j(y^-x^yJ 

Cx  =  j(y^-Z^) 

Cy  =  J  (x-qZ^  ~  x^zti) 

Cz  =  j(x^yq-y^) 
x-t^x  ~  y-t^y  “  Z^z 
Tit  =  ~  XtTlx  —  y^Tly  —  Z^Tlx 
Ct  ~  ~  XtCx  “  y-tCy  “  ZtCz 


where  J  is  the  Jacobian  of  the  transformation  and  is  given  by 


J=  1/J-1 


j-i  _  d(x,y,z,t)  _ 

afen.i;,*) 


10  0  0 
X-t  X^  Xyj  Xj^ 

yxy^yqy; 

Zt  Z^  z^  z^ 


j-1  =  x^z^  +  x^y^  +  x^y^  -  x^y^  -  x^y^  -  x^z^ 


(A4) 


(A5) 


(A6) 


(A7) 
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If  the  transformation  is  a  result  of  grid  generation  then  the  metrics  and  the  Jacobian  can  be 
computed  numerically  using  central  differences. 

Vinokur  (48)  showed  that  the  governing  equations  can  be  transformed  into  strong 
conservation  law  form,  by  first  dividing  Equation  A3  by  the  Jacobian  and  then  using  the 
chain  rule  to  bring  the  Jacobian  inside  the  differential  operators.  Equation  A3  becomes 


the  last  four  terms  on  the  left  hand  side  and  the  last  three  terms  of  the  right  hand  side  are 
known  as  the  invarients  of  the  transformation.  It  can  be  shown,  by  using  the  definition  of 
the  metrics,  that  these  terms  are  zero  (32:6). 

The  contravarient  velocities  U,  V,  and  W  are  defined  to  be  the  velocities  in 
directions  normal  to  planes  of  constant  Tj,  and  t,  respectively  and  are  given  by 

U  +  +^yV  +  ^zW 

V  =  T|t  +  T)xU  +  TJyV  +  T|ZW  (A9) 

w  =  Ct  +  Cxti  +  Cy  V  + 
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Applying  these  definitions  to  Equation  A8,  Pullian  (38:159-160)  showed  that  the  governing 
equations  written  in  strong  conservation-law  form  and  generalized  coordinates  are  given  by 


8tQ  +  dsE  +  9-r]  F  +  d^G  =  (dj;  Ev  +  d-^  v  +  d^Gv) 


(A10) 


Q  =  Jl  pv 


(All) 


pu  i  r  pv  '  r  Pw  i 

puU  +  ^xp  puV  +  r|xp  puW  +  CxP 

E  =  J_1  pvU  +  ^yP  F  =  Jl  pvV  +  qyp  G=J"'  pvW  +  £yp  (A12) 

pwU  +  £zp  pwV  +  qzp  pwW  +  £zp 


.(e+p)U+^tp 


Ev  =  J~l 


Fv  =  J-1 


_(e+p)V  +  ritp 


^x^xx  +  ^y^xy  +  ^z'txz 
^x^yx  +  ^y'tyy  +  ^z^yz 

^x^zx  +  4y^zy  +  4z^zz 
4xPx  “*■  ^yPy  +  ^zPz 

0 

TlxXxx  +  %Xxy  +  T|zXxz 
Tlx’Cyx  +  Hy^yy  +  Hz^yz 
Tlx*zx  +  TlyXzy  +  Tlzt  zz 
rixPx  +  T)y^y  +  qzpz 


(e+p)W  +  ClP_ 


(A  13) 
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Gv  =  J-1 


Cx^xx  +  Cy^xy  +  Cz^xz 
Cx^yx  +  Cy^yy  +  Cz^yz 


Cx^zx  +  ^yl-zy  +  Cz^zz 
L  CxPx  +  CyPy  +  CzPz  _ 


Px  =  Y(iPr_1  3xei  +  ux„  +  vxxy  +  wtxz 

Py  =  Y^lPr-1  dyC[  +  Ulyn  +  V  Tyy  +  WTyZ 

P2  =  YM-Pr-1  Ozei  +  utzx  +  vxzy  +  wx^ 


and  the  cartesian  derivatives  become 


«x  =  5xU$  +  Tlxll,,  +  CxU; 

Uy  =  £yU$  +  TJyU^  +  £yU? 
UZ  =  £ZU^  +  +  CzU; 

Vx  =  ^  +  r\x\  +  CxVj; 

Vy  =  %yV%  +  HyV^  +  CyV^ 

Vz  =  +  T\zVf|  +  CzV; 

Wx  =  CXW^  +  TlxW^  +  Cx  W(; 

Wy  =  ^yW^  +  TlyW^  +  ^yW^ 
Wz  =  ^zW^  +  T]z  +  CzW^ 


(A14) 


(A  15) 
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Pulliam  has  shown  (  32:99-100)  that  the  invisicid  flux  Jacobians,  A,  B,  and  C,  are 
obtained  by  time  linearizations  of  the  inviscid  flux  vectors,  E,  F,  and  G,  and  are  given  by 


A 


3E 

3Q 


^-uQi 

5t  +  01  -UY-2)u 

£,y())2  -  V01 

^xV-4y(Y-  l)u 

^z4>2-we! 

^xW-Uy-  1)u 

-9i(ye/p-202) 

Sx¥~(Y-  1)u«i 

^z 

0 

^yU-^x(Y-l)v 

^zu-cx(y-  l)w 

Uy-i) 

^,  +  0a-^y(Y-2)v 

^zV-^y(Y-l)w 

5y(Y“l) 

^yW-UY-  l)v 

^t  +  0a-^z(Y-2)w 

Uy-i) 

^y¥-(Y-  1)V01 

^z¥-{Y-l)w6i 

^t  +  Y0i 

Til 

flx 

Tlx<t>2  -  U02b 

flt  +  02-flx(Y-2)u 

Tly<t>2  -  V02 

flxV-T1y(Y-  l)u 

flz<l>2  -  W02 

flxW-Tlz(Y-  l)u 

-02(Ye/p  -  2p2) 

flx¥-(Y-  l)u02 

fly 

flz 

0 

flyU-Tlx(Y-  l)v 

flzu-flx(Y-  l)w 

flx(Y-  1) 

Ti,  +  02-fly(Y-2)v  TizV-riy(Y-l)w 

fly(Y-l) 

flyW-Tlz(Y-  l)v 

flt  +  02-flz(Y-2)w 

flz(Y-  1) 

fly¥-(Y-  l)v02 

flz¥-(Y-  1)W02 

flt  +  Y02 

(Bl) 


(B2) 


dQ 


e« 

Cx 

C-/-U03 

Ct  +  03-Cx(Y-2)u 

Cy4>2-v03 

CxV-Cy(Y-  Ou 

^  -  W03 

CxW-Cz(Y-  1)u 

-03(ye/p  -  2<J>2) 

^xV-(Y-  l)u03 

Cz 

0 

CyU-Cx(Y-  l)v 

^zU-Cx(Y~  ])w 

Cx(Y-D 

Ct  +  03-Cy(Y-2)v 

CzV-CyfY-Uw 

Cy(Y-l) 

CyW-Cz(Y-  l)v 

Ct  +  03-Cz(Y-2)w 

Cz(Y-l) 

CyV-(Y"  1)V03 

CzV  —  (Y“  1  )w03 

Ct  +  Y03 

(B3) 


where 


(B4) 


V  =  -ye/p  -  <J> 


(B5) 


9l  =  ^XU  +  ^yV  + 

02  =  +  %v  +  ruw 

03  *  CxU  +  Cyv  +  CzW 


(B6) 


The  viscous  flux  Jacobian  (M)  is  obtained  from  the  linearization  of  the  viscous  flux 
vector  (S).  Using  a  Taylor  series,  Steger  has  shown  the  Jacobian  is  given  by  (43:3-4) 


bO 


0  0  0  0  0 

m2i  atS^p-1)  a25^p_1)  a^p-1)  0 

M  =  ^| -  =  m31  a28c(p_1)  cuS^p"1)  asS^p-1)  0  (B7) 

3Q 

m4i  a-3S^(p- 1 )  a55^p_1)  a65^p-!)  0 

msi  m52  msi  ms4  ao8^(p~ 1 )  _ 

m2i  =  ai8^{-u/p)  +  a2S^-v/p)  +  a35^(-w/p) 
m3i  =  a28^-u/p)  +  a4S^-v/p)  +  as5^(-w/p) 
nui  =  a36^(-u/p)  +  a58^-v/p)  +  a65^(-w/p) 
msi  =  ai6^{-u2/p)  +  a25;;(-2uv/p)  +  a38^-2uw/p)  +  a48^(-v2/p)  + 

<*68^-w2/p)  +  oco8(j[-e/p2)  +  aoS<{(u2  +  v2  +  w2Vp]  (B8) 

m52  =  -m2i  -  aoS^u/p) 
m53  =  -m3i  -  ctoS^v/p) 

m54  =  -nxu  -  a08$(w/p) 

cx0  =  7pPr1(;5  +  Cy  +  Cz) 
a,  =  p[(4/3)d  +  C?  +  d] 
a2  =  (p/3)CxCy 

a2  =  (p/3)CxCz  (B9) 

«4  =  p[d  +  (4/3)Cy  +  d] 
a5  =  (M/3)CyCz 
a6»n[d  +  ^  +  (4/3)d] 
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Warming,  Beam,  and  Hyett  (44:1037-1041)  have  shown  that  the  invisicid  flux 
Jacobians  can  be  diagonalized  since  they  have  real  eigenvalues  and  a  complete  set  of 
eigenvectors.  The  flux  Jacobians  are  written  in  terms  of  their  eigenvalues  and  eigenvectors 
as 


-l\n 


A=|T;A',T. 

B"  =(TnA„T-')" 

P-kf 


(Cl) 


where  is  the  eigenvector  matrix  of  A  and  likewise  for  B  and  Tt;  for  C.  The 


eigenvalue  matrices  are  given  by 

r 

u 

0 

0 

0 

0  ‘ 

0 

u 

0 

0 

0 

0 

0 

u 

0 

0 

0 

0 

0 

U+aKj 

0 

- 

0 

0 

0 

0 

U-aici  _ 

- 

V 

0 

0 

0 

0  1 

0 

V 

0 

0 

0 

> 

ii 

0 

0 

V 

0 

0 

0 

0 

0 

V+aK2 

0 

- 

0 

0 

0 

0 

V-aic2  J 

- 

w 

0 

0 

0 

0  “ 

0 

w 

0 

0 

0 

A?- 

0 

0 

w 

0 

0 

0 

0 

0 

W+a*3 

0 

0 

0 

0 

0 

W-aK3 

where  a  is  a  characteristic  speed  and 


K2  =  ('li  +  'l?  +  ril)1'2 

K3=(d  +  C5  +  S2)1/2 


The  eigenvector  matrices  are 


jx 

Ts= 

ixV+^zp 

^yu-4zP 
~  ^yVL 

SxW-^yP 

_  kx02/(Y-l)  +  p(iv-qyw) 

kyW+lxP 
Jy<t>2/(Y-1)  +  p(£>P 

*-£zu). 

a 

a 

izU+5yp 

a(u+£xa) 

a(u-^xa) 

&-lxP 

a(v+|yaj 

a(v-?ya) 

*** 

a(w+£za) 

a(w-£za) 

Ez«l>2/(Y-l)  +  p(^yU-^v)]  a[(«j)2 

-r  a2)/(y-l)  +  0ia] 

a[(<(>2  +  a2)/(y-l) 

T,= 


Tlx 

Tly 

TlxU 

TlyU-rizP 

TlxV+Tlzp 

rfyv 

TlxW-TlyP 

%W  +  T1xP 

L  [t1x02/(Y-  1 )  +  p(riuv-riyw)]  [rj|’y<J)2/(y—  1 )  +  p(rfxW-T)zu)] 


Tlz 

a 

a 

TlzU+TlyP 

a(u+r[xa) 

a(u-rfxa) 

TlzV-TlxP 

a(v+nya) 

afv-iiya) 

TlzW 

a(w+riza) 

a(w-r[za) 

r-l)  +  pfrfyU-Tlxv)] 

a[(<j>2  +  a2)/(Y-l)  +  02a 

a[(<|)2  +  a2)/(Y-l)-02a 

L  [u2A7-D  +  p(CzV-Cyw)]  ky02/(Y-l)  +  p(CxW-Czu) 

\  °L  a~ 

CzU+CyP  a(u+C*a)  a(u-Cxa) 

CzV-CxP  a(v+Jya)  a(v-Jya) 

Czw  a(w+£za)  a(w-<za) 

l^z<!>2/(Y-l)  +  p(CyU-Cxv)]  ajjp2  +  a2)/(y—  1 )  +  83a]  a[(<J>2  +  a2)/(y-l)-03a 


5x(  1  -02/a2)  -  fezv-qy  w)/p 

^x(Y-l)u/a2 

£y(  1  -4>2/a2)  -  (^XW -^zu)/p 

Jy(Y-l)u/a2+Jx/p 

Ui-4>2/a2)-kyU-qxv^p 

4z(Y-I)u/a2+Vp 

p(<l>2-0ia) 

-^Y-^-^xa 

(3U2+0ia) 

-WY-l^+^a 

^x(Y-l)v/a2+^p 

^x(Y-l)w/a2~^y/p 

-UY-iya2 

^y(Y-l)v/a2 

^Y-l)w/a2+^x/p 

-^Y-iya2 

UY-1  )v/a2+£x/p 

UY_1Wa2 

-UY-^a2 

-P[(Y~  1  K-^ya] 

-p[(Y-l)w-cza] 

P(Y-D 

-p[(Y-l)v+$ya] 

-P[(Y-l)w+^za] 

P(Y-l) 

rfx(  1  — 4>2/a2)  -  Rzv-riywyp 

rfx(Y-l)u/a2 

%( 1  -<>2/a2)  -  (ruw-nzuyp 

Tly(Y— 1  )u/a2  -t-  tTx/P 

T\zi  1  -02/a2)  -  fiyU-Tlxvj/p 

rf^Y-l^+VP 

p(<|>2-02a) 

-(3(Y-l)u-Tlxa 

p(<j>2+02a) 

-P(Y-l)u+rtxa 

The  use  of  the  solid  state  disk  allows  for  the  calculations  to  be  performed  on  grids 
that  require  more  computer  memory  than  is  available.  The  purpose  of  this  appendix  is  to 
evaluate  the  overhead  associated  with  using  the  solid  state  disk.  The  code  was  run  for  50 
iterations  in  order  to  get  a  representative  feel  for  differences  in  central  processor  unit  (CPU) 
times,  input/output  (I/O)  requirements,  and  job  size. 

A  grid  was  generated  for  a  delta  wing  that  was  small  enough  to  fit  entirely  in  the 
core  of  the  Cray  XMP-12  computer.  The  grid  consisted  of  32,000  points  which  is  about 
the  maximum  number  of  points  that  can  be  processed  at  one  time.  The  grid  had  20  points 
in  the  streamwise  direction,  40  points  in  the  spanwise  direction,  and  40  points  in  the 
normal  direction. 

The  accuracy  of  the  solution  is  not  affected  by  the  choice  of  methods.  The  only 
difference  in  the  operation  of  the  code  is,  the  solid  state  disk  loads  and  unloads  planes  of 
data  as  needed  for  calculations.  This  resulted  in  a  40  percent  increase  in  CPU  time.  This 
number  will  vary  depending  on  the  number  of  data  planes  that  are  processed  per  iteration. 

The  in-core  run  required  a  job  size  of  1 ,420,000  words  to  execute.  This  used 
approximately  80  percent  of  the  available  core  memory.  The  SSD  run  required  260,000 
words  of  core  to  execute.  The  problem  encountered  when  running  a  job  that  requires  a 
large  block  of  core  memory  is  that  the  job  may  sit  for  extended  periods.  The  in-core  job  sat 
for  almost  30  hours,  while  the  out-of-core  job  sat  for  about  5  minutes,  and  was  finished  in 
about  30  minutes. 

The  last  area  examined  was  the  cost  of  running  the  job.  The  in-core  cost  was  $40 
dollars,  while  the  out-of-core  jobcost  was  $375.  It  is  believed  that  this  cost  is  due  to  an 
input/output  charge.  The  SSD  made  almost  250,000  I/O  requests,  while  the  in-core  had 
150.  If  this  is  the  case,  then  using  the  SSD  may  become  prohibitive  from  a  financial  point 
of  view. 
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This  appendix  contains  the  job  control  language  and  parameter  input  files  necessary 


to  run  the  ARC3D  code  on  the  Cray  XMP-12  computer  at  Wright-Patterson  AFB,Ohio. 

The  first  file  is  used  for  compiling  the  code  and  storing  it  on  the  Cray.  The  second  file  is 
the  file  used  for  submitting  a  job  to  the  Cray.  The  code  will  look  for  three  files  before  it 
begins  execution,  1)  a  restart  file,  if  requested,  2)  a  grid  file,  and  3)  a  parameter  input  file. 
The  third  file  is  an  example  of  a  parameter  input  file.  Upon  termination  the  program  returns 
a  restart  file  and  a  parameter  output  file.  Due  to  large  amount  of  CPU  time  required,  the 
code  was  run  in  batches  of  200  time  iterations. 


*********************  'J'q  Compile  ARC3D  Code  ************************* 

JOB,JN=CRANKCMP,T=50,MFL= 1 50000.  SMITH- AFIT/Em  -54731 

ACCOUNT,  AC=xxxxxxx,APW=xxxxxx,US=xxxxxxx,UPW=xxxxxx. 

* 

THIS  FILE  IS  USED  TO  COMPILE  THE  ARC3D  CODE  ON  THE  CRAY  XMP 

*.  A  COPY  OF  THE  COMPILED  PROGRAM  IS  STORED  ON  THE  CRAY. 

* 

PURGE  THE  PREVIOUS  COPY. 

* 

ACCESS,DN=ARC3D,ID=D880165,UQ. 

DELETE, DN=ARC3D. 

RELEASE,DN=ARC3D. 

* 

*.  GET  SOURCE  CODE  FROM  THE  CFS. 

*_ 

FETCH,DN=CODE,SDN= ARC3  D  ,MF=CB  ,DF =CB  ,A 

TEXT=’SREAD,ARC3D,ARC3D.CTASK,ALL.’. 

* 

*.  COMPILE  CODE  AND  WRITE  TO  BINCODE 
* 

CFT,I=CODE,B=BINCODE,ON=Z,OPT=FULLIFCON,MAXBLOCK=4000,L=0. 

* 

*  *  SAVE  EXECUTABLE  ON  THE  CRAY 

* 

S  A  VE,DN=B  INCODE, PDN=ARC3D,ID=D880 1 65. 

* 

SEND  BACKUP  TO  THE  CFS 

* 

*.  DISPOSE,DN=BINCODE,SDN=ARC3DEXE,MF=CB,DC=ST,DF=TR,A 
*.  TEXT=’CTASKrALL.SWRITE,ARC3DEXE,ARC3DEXE'. 

/EOF 


I 


*************************  ARC3D  Input  File  **************************** 

JOB  JN=CRANK,T=  1 600,MFL=500000,SSD=  1 1000,CL=P2.  SMITH/54731 

ACCOUNT, AC=xxxxxxx,APW=xxxxxx,US=xxxxxxx,UPW=xxxxxx. 

* 

THIS  FILE  RUNS  THE  ARC3D  COMPUTER  CODE  ON  THE  ASD 

*.  CRAY  XMP  USING  THE  SOLID  STATE  DISK. 

* 

*  SET  UP  FILES  ON  THE  SSD. 

* 

ASSIGN,DN=FT1 1  ,D  V=SSD-0-20,RDM,U. 

ASSIGN,DN=FT12,DV=SSD-0-20,RDM,U. 

* 

*.  GET  THE  RESTART  FILE  FROM  THE  CFS. 

* 

FETCH,DN=RESTIN,SDN=RESTIN,MF=CB,DF=TR,A 

TEXT='SREAD,RESTIN,RESTIN.CTASK,ALL.'. 

* 

GET  THE  GRID  FROM  THE  CFS. 

* 

FETCH,DN==GRIDIN,SDN=BNGRID,MF=CB,DF=TR,A 

TEXT='SREAD,BNGRID,BINGRID.CTASK,ALL.,. 

* 

*.  GET  INPUT  FILE  FROM  THE  CYBER. 

*_ 

FETCH,DN=PARAIN,SDN=PARAIN,MF=CB,DF=CB,A 

TEXT='GET,PARAIN.CTASK,ALL.'. 

* 

*.  GET  COMPILED  ARC3D  CODE  FROM  THE  CRAY. 

*_ 

ACCESS,DN=$BLD,PDN=ARC3D,ID=D880165. 

OPTION,STAT=ON. 

* 

LOAD  AND  RUN. 

* 

LDR. 

*  _ 

SAVE  THE  RESTART  TO  THE  CFS. 

* 

DISPOSE,DN=RESTOUT,SDN=RESTOUT,MF=CB,DC=ST,DF=TR,DEFER,A 
TEXT=’CTASK,ALL.SWRITE, RESTOUT, RESTOUT. 

SAVE  THE  GRID  OUTPUT  TO  THE  CFS. 

* 

*.  DISPOSE,DN=GRIDOUT,SDN=GRIDPLOT,MF=CB,DC=ST,DF=TR,A 
*.  TEXT=’CTASK,ALL.SWRITE,GRIDPLOT,GRIDPLOT. 

* 

*.  BATCH  PARAMETER  OUTPUT  FILE  TO  THE  WAIT  QUEUE. 

* 

DISPOSE,DN=PAROUT,SDN=PAROUT,MF=CB,DC=ST,DF=CB,DEFER,A 

TEXT=’CTASK,ALL.ROUTE,PAROUT,DC=WT,UJN=PAROUT,UN=D880165.’. 
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EXIT. 

DISPOSE,DN=RESTOUT,SDN=RESTOUT,MF=CB,DC=ST,DF=TR,DEFER  ,A 

TEXT='CTASK,ALL.SWRITE,RESTOUT,RESTOUT’. 

DISPOSE,DN=PAROUT,SDN=PAROUT,MF=CB,DC=ST,DF=CB,DEFER,A 

TEXT='CTASK,ALL.ROUTE,PAROUT,DC=WT,UJN=PAROUT,UN=D880165.'. 

DUMPJOB. 

DEBUG. 

* 

/EOF 


*********************  Parameter  Input  File  ************************ 


F 

28,65,70 

TT 

TTF 

50 

T 

1 

21  25 
T 

2 

24  25 
24  24 
T 
1 

24  21 

TF 

TF 

686650,  0.72,  234.4 
2.5,10.,0. 

1 

3 

0.50,0.75 

0.50,0.75 

0.50,0.75 

1 

F 

4 


41  25 
1.0 
2 
1 


3  29 


2.00  200 


INCORE  (T/F) 

JMAX,KMAX,LMAX 
RESTART.STORE 
READGRID, READ  ALL,  WRTTGRID 
NP 

PJLINE  (ENTER  NJLINE(<1 1),THEN  EACH  PAIR  OF  K,L) 
NJLINE 

JLINK(NJ),JLINL(NJ) 

PKLINE  (ENTER  NKLINE(<  1 1 )  .THEN  EACH  PAIR  OF  J,L) 
NKLINE 

KLIN  J(NK)  ,KLINL(  NK) 

PLLINE  (ENTER  NLLINE(<1 1),THEN  EACH  PAIR  OF  J,K) 
NLLINE 

LLINJ(NL),LLINK(NL) 

IFSCOR,METAV 

VISCOUS, TURBULNT  (IF  VISCOUS  THEN  ENTER 
RE,PR,TINF) 

RE,PR,TINF 
FSMACH,  ALP,  Y  AW 
METH 

ISPEC  ARTIFICIAL  VISCOSITY  STUFF 

DIS2X.DIS4X 

DIS2Y.DIS4Y 

DIS2Z  DIS4Z 

IVARDT  0-CONST  1 -VARIABLE  TIME  STEP 
PERIODIC  IN  K 

IBC(IF  IBC=3  THEN  ENTER  JTAIL1,JTAIL2) 

(THERE  IS  NO  4!!) 

KEDGE  LSURF  JSTART  JEND  IBC=4  STUFF 

CLENGTH  (CHECK  THIS  NUMBER) 

IORDER 

NUMDT  (NUMBER  OF  TIMESTEP  SEQUENCES  TO 
FOLLOW) 

DTSEQ(I),ITERDT(I)  (TIME  STEP  AND  NUMBER  OF 
TIMESTEPS) 


/EOF 


This  program  was  originally  written  by  Dr.  Phil  Webster  of  the  Air  Force  Flight 
Dynamics  Laboratory  to  generate  a  grid  for  a  delta  wing.  It  has  since  been  modified  by  the 
author  to  generate  the  grid  for  a  cranked  delta  wing.  A  branch  cut  between  the  upper  and 
lower  surfaces  of  the  body  allows  for  the  use  of  an  H  grid.  This  approach  was  used  to 
eliminate  the  small  radius  of  curvature  generated  at  the  leading  edge  using  a  C  grid.  Points 
are  clustered  in  high  gradient  areas  using  an  exponential  stretching  function.  There  is  no 
provision  for  a  wake  region,  therefore  this  code  is  only  used  to  generate  grids  for 
supersonic  flows. 


PROGRAM  DELGRD 
C 

C  THIS  PROGRAM  WAS  MODIFIED  TO  CLUSTER  POINTS  AT  THE  WING 
C  LEADING  EDGE  BY  USING  ELLIPSES  FOR  THE  UPPER  AND  LOWER 
C  BOUNDARIES  OF  THE  DOMAIN.  FRS  8  AUG  88 
C 

C  THIS  PROGRAM  WAS  MODIFIED  TO  GENERATE  A  GRID  FOR  A  CRANKED 
C  DELTA  WING  CONFIGURATION.  FRS  28  AUG  88 
C 

C  PROGRAM  GENERATES  A  GRID  (HOPEFULLY)  TO  FIT  THE  DELTA  WING 
C  MODIFIED  FOR  A  DELTA  WING  WITH  NO  WAKE  REGION 
C 

CDIR$  NOVECTOR 
REALLEN 

COMMON  X(28,65,70),Y(28, 65,70), Z(28, 65,70) 

LOGICAL  LSTOP,LCHECK,LCRANK 
LCHECK  =  .TRUE. 

LSTOP  =  .FALSE. 

LCRANK  =  .TRUE. 

C 

C  OUTPUT  FILES  FOR  IRIS  (MUST  BE  COMMENTED  OUT) 

C 

C  OPEN  (UNIT=10,FILE='GRID.BIN\FORM=BINARY’) 

C  OPEN  (UNTT=9,FILE='GRID.DAT) 

C 

C  FILES  FOR  THE  KIRTLAND  AFB  CRAY 
C 

OPEN  (10,FTLE=’CGRID',STATUS=’NEW',FORM='UNFORMATTED') 

IF  (LCHECK)  OPEN  (9,FILE=’CKGRID’,STATUS='NEW’, 

&  FORM=’FORMATTED') 

PI  =  ASIN(l.O)  *  2.0 
LEN  =  .100 


nnn  o  non 


ALPHA  =  10.0 

ALPHA  =  (ALPHA  /  180.0)  *  PI 
OMEGA  =40.0 

OMEGA  =  (OMEGA  /  180.0)  *  PI 
DELTA  =  1.0E-3  *  LEN 
XSTART  =  0.05  *  LEN 
XEND  =  XSTART  +  LEN 
XTOTAL  =  XEND 
JSTART  =  3 
KMID  =  4 


INPUT  GEOMETRY  FOR  DELTA  WING 

IF(.NOT.LCRANK)THEN 
SWEEP  =  80.0 

SWEEP  =  ((90.0  -  SWEEP)  /  180.0)  *  PI 
BASE  =  LEN  *  TAN(SWEEP) 

XMID  =  XSTART  +  0.90  *  LEN 


NJ  =  30 

JMID  =  28 

JEND  =  30 

NK=55 

KEDGE=34 

NL=65 

LSURF=21 

INPUT  GEOMETRY  FOR  CRANKED  DELTA  WING 
ELSE 

XCRANK=.402*LEN+XST  ART 
SWEEP  1=80 

SWEEP1=((90.-SWEEP1)/180.)*PI 

SWEEP2=69 

SWEEP2=((90.-SWEEP2)/180.)*PI 

BASE  1 =(XCRANK~XSTART)*T  AN(S' WEEP  1 ) 

BASE2=(XEND-XCRANK)*TAN(SWEEP2)+BASE1 

XMID1=(XCRANK-XSTART)*.70+XSTART 

XMID2=XSTART+.70*LEN 

NJ=28 
JMID  1=8 
JCRANK=13 
JMID2=24 
JEND=28 
NK=65 
KEDGE=41 
NL=70 
LSURF=25 
END  IF 

IF  (LCHECK)  THEN 

WRITE  (9,*)  ’  DELTA  WING  LENGTH  IS  \LEN 


7/ 


IF(.NOT.LCRANK)THEN 

WRITE  (9,*) '  HALF  THE  WING  ROOT  CHORD  IS  ’.BASE 
THICK=BASE*TAN(OMEGA) 

WRITE  (9,*) '  ROOT  THICKNESS  IS  V THICK 
ELSE 

WRITE  (9,*) '  HALF  THE  WING  ROOT  CHORD  IS  \BASE2 
THICK=B  ASE2*T  AN(OMEG  A) 

WRITE  (9,*) '  ROOT  THICKNESS  IS  '.THICK 
END  IF 
ENDIF 
C 

DO  1  J  =  1,NJ 
DO  1  K  =  1.NK 
DO  1  L  =  1.NL 
X(J,K,L)  *  0.0 
Y(J,K,L)  =  0.0 
Z(J,K,L)  =  0.0 
1  CONTINUE 
C 

C  CONSTANTS  FOR  THE  CALCULATION  OF  THE  LIMITS  OF  THE  DOMAIN 
C 

YMAX1  =0.15  *  LEN 
ZMAX1T  =  0.15  *  LEN 
ZMAX1B  = -0.15  *  LEN 
IF(.NOT.LCRANK)THEN 
YMAX2  =  0.85*  LEN 
ZMAX2T  =  1.05*  LEN 
ZMAX2B  =  -0.65  *  LEN 
ELSE 

YMAX2  =  0.95*  LEN 
ZMAX2T  =  1.05  *  LEN 
ZMAX2B  =  -0.80  *  LEN 
END  IF 

XTOTSQ  =  XTOTAL  *  XTOTAL 
XTOTCU  =  XTOTSQ  *  XTOTAL 
C 

C  THE  CONSTANT  CURVE  CONTROLS  THE  AMOUNT  OF  CURVATURE  IN  THE 
C  MAX  AND  MIN  Z  PLANES.  CURVE  MUST  BE  GREATER  THAN  OR  EQUAL 
C  TO  0  AND  LESS  THAN  1 .  IF  CURVE  IS  EQUAL  TO  0  A  RECTANGULAR 
C  GRID  RESULTS.  IF  CURVE  EQUALS  1  A  SINGULARITY  OCCURS  AT 
C  Y=YMAX. 

C 

CURVE=0. 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  X  DIRECTION 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  DIST  X  POINTS  ALONG  THE  L  =  LSURF  AND  K  =  l 
C 

C  AHEAD  OF  BODY 
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IF  (LCHECK)  WRITE  (9,*) '  X  DIRECTION' 
J  =  1 

DO  90  K  =  1,NK 
DO  90  L  =  1,NL 
X(J,K,L)  =  O.o 
90  CONTINUE 

IF  (LCHECK)  WRITE  (9,*)  J,X(J,1,LSURF) 
C 

NUMB  =  JSTART  -  1 

DIST  =  XSTART 

DETA  =1.0/  FLOAT(NUMB) 

DX  =  20.0*DELTA 

CALL  FINDC(DX,DIST,NUMB  ,C,  1 ) 

DO  100  J  =  2, JSTART- 1 
JJ  =  JSTART  -  J 
ETA  =  DETA  *  FLOAT(JJ) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1 .0) 
XTEMP  =  XSTART  -  FAC  *  DIST 
DO  99  K  =  1,NK 
DO  99  L  =  1,NL 
X(J,K,L)  =  XTEMP 

99  CONTINUE 

IF  (LCHECK)  THEN 
Q  =  (XSTART  -  XTEMP)  /  DX 
WRITE  (9,*)  J,X( J,  1  ,LS URF) ,Q 
ENDIF 

100  CONTINUE 
C 

J  =  JSTART 
DO  101  K=  1,NK 
DO  101  L=  1,NL 
X(J,K,L)  =  XSTART 

101  CONTINUE 

IF  (LCHECK)  WRITE  (9,*)  J,X(J,1,LSURF) 

IF(.NOT.LCRANK)THEN 

FIRST  SECTION  OF  BODY 

NUMB  =  JMID  -  JSTART 
DIST  =  XMID  -  XSTART 
DETA  =1.0/  FLOAT(NUMB) 

DX  =  15.0  *  DELTA 
CALL  FENDC(DX,DIST,NUMB ,C,  1 ) 

DO  1 10  J  =  JSTART+1  ,JMID- 1 
JJ  =  J  -  JSTART 
ETA  =  DETA  *  FLOAT(JJ) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1.0) 
XTEMP  =  FAC  *  DIST  +  XSTART 
DO  109  k  =  1,NK 
DO  109  L=  1,NL 
X(J,K,L)  =  XTEMP 
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109  CONTINUE 

IF  (LCHECK)  THEN 
Q  =  (XTEMP  -  XSTART)  /  DX 
WRITE  (9,*)  J,X(J,1,LSURF),Q 
ENDIF 

1 10  CONTINUE 
C 

J  =  JMID 
DO  111  K=  1,NK 
DO  111  L=  1,NL 
X(J,K,L)  =  XMID 

111  CONTINUE 

IF  (LCHECK)  WRITE  (9,*)  J,X(J,1,LSURF) 
NOW  FOR  THE  SECOND  PART  OF  THE  BODY 


NUMB  =  JEND  -  JMID 
DIST  =  XEND  -  XMID 
DETA  =1.0/  FLO  AT (NUMB) 

DX  =  5.0  *  DELTA 

CALL  FINDC(DX,DIST,NUMB,C,1) 

DO  120  J  =  JMID+1 JEND-1 
JJ  =  JEND  -  J 
ETA  =  DETA  *  FLOAT(JJ) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1 .0) 
XTEMP  =  XEND  -  FAC  *  DIST 
DO  119  K=  1,NK 
DO  1 19  L  =  1,NL 
X(J,K,L)  =  XTEMP 

19  CONTINUE 

IF  (LCHECK)  THEN 
Q  =  (XEND  -  XTEMP)  /  DX 
WRITE  (9,*)  J,X(J,1,LSURF),Q 
ENDIF 

20  CONTINUE 


J  =  JEND 
DO  121  K=  1,NK 
DO  121  L=  1,NL 
X(J,K,L)  =  XEND 
[21  CONTINUE 

IF  (LCHECK)  WRITE  (9*)  J,X(J,1,LSURF) 

DO  THE  CRANKED  WING 
ELSE 

DISTRIBUTE  POINTS  FROM  JSTART  TO  JMID1 

NUMB  =  JMID  1  -JSTART 
DIST  =  XMJD1  -  XSTART 
DETA  =1.0/  FLOAT(NUMB) 

DX  =  50.0  *  DELTA 


non  n  non 


CALL  FINDC(DX,DIST,NUMB,C,  1 ) 

DO  125  J  =  JSTART+1 JMID1-1 
JJ  =  J  -  JSTART 
ETA  =  DETA  *  FLOAT(JJ) 

FAC  =  (EXP(C*ETA)  -  1.0)/  (EXP(C)  -  1.0) 
XTEMP  =  FAC  *  DIST  +  XSTART 
DO  124  K  =  1,NK 
DO  124  L=  1,NL 

124  X(J,K,L)  =  XTEMP 
IF  (LCHECK)  THEN 

Q  =  (XTEMP  -  XSTART)  /  DX 
WRITE  (9,*)  J,X(J,1,LSURF),Q 
END  IF 

125  CONTINUE 
J  =  JMID1 

DO  126  K=  1,NK 
DO  126  L  =  1,NL 

126  X(J,K,L)  -  XMID 1 

IF  (LCHECK)  WRITE  (9,*)  J,X(J,1,LSURF) 

DISTRIBUTE  POINTS  FROM  JMID1  TO  J CRANK 

NUMB  =  JCRANK-JMID1 
DETA  =1.0/  FLOAT  (NUMB ) 

DX  =  8.0  *  DELTA 

DIST  =  XCRANK-  .5  *DX-XMID  1 

CALL  FINDC(DX,DIST,NUMB,C,  1 ) 

DO  130  J  =  JMID 1+1 ,  JCRANK- 1 
JJ  =  JCRANK  -  J 
ETA  =  DETA  *  FLOAT(JJ) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1 .0) 
XTEMP  =  XCRANK-.5*DX-FAC*DIST 
DO  129  K  =  1,NK 
DO  129  L  =  1,NL 

129  X(J,K,L)  =  XTEMP 
IF  (LCHECK)  THEN 

Q  =  (XCRANK-.5*DX-XTEMP)/DX 
WRITE  (9,*)  J,X(J,1,LSURF),Q 
END  IF 

130  CONTINUE 

J  =  JCRANK 
DO  131  K  =  1,NK 
DO  131  L=  l,NL 
X(J,K,L)  =  XCRANK-. 5*DX 

131  X(J+1,K,L)  =  XCRANK+.5*DX 
IF  (LCHECK)  THEN 

WRITE(9,*)J,X(J,1,LSURF) 

WRITE(9,*)J+1  vX(J-i- 1 , 1  ,LSURF) 

END  IF 

DISTRIBUTE  POINTS  FROM  JCRANK+1  TO  JMID2 
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NUMB  =  JMID2  -  (JCRANK+1) 

DIST  =  XMID2  -  X( JCRANK+ 1 , 1  ,LSURF) 

DETA  =1.0/  FLOAT(NUMB) 

DX  =  8.0  *  DELTA 

CALL  FINDCl  DX,DIST,NUMB  ,C,  1 ) 

DO  135  J  =  JCRANK+2JMID2-1 
JJ  =  J  -  (JCRANK+1) 

ETA  =  DETA  *  FLOAT(JJ) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1.0) 
XTEMP  =  FAC  *  DIST  +  X(JCRANK+ 1,1, LSURF) 
DO  136  K=  1,NK 
DO  136  L=  1,NL 

136  X(J,K,L)  =  XTEMP 
IF  (LCHECK)  THEN 

Q  =  (XTEMP  -  X(JCRANK+ 1,1, LSURF)  )/DX 
WRITE  (9,*)  J,X(J,  1  ,LS URF) ,Q 
ENDIF 

135  CONTINUE 

J  =  JMID2 
DO  137  K  =  1,NK 
DO  137  L  =  1,NL 

137  X(J,K,L)  =  XMID2 

IF  (LCHECK)  WRITE  (9,*)  J,X(J,1,LSURF) 

DISTRIBUTE  POINTS  FROM  JMID2  TO  JEND 

DX  =  5.0*  DELTA 
DO  140  K=  1,NK 
DO  140  L  =  1,NL 
X(J+1,K,L)  =  XMID2+0.08  *  LEN 
X(J+2,K,L)  =  XMID2+0. 18*LEN 

140  X(J+3,K,L)  =  XMID2+0.30*LEN-DX 
IF(LCHECK)THEN 

DO  141  J =JMID2+ 1  JEND- 1 
Q=(XEND-X(J,  1  ,LSURF))/DX 
WRITE(9,*)J,X(J,1,LSURF),Q 

141  CONTINUE 
ENDIF 

J  =  JEND 
DO  142  K  =  1,NK 
DO  142  L  =  1,NL 

142  X(J,K,L)  =  XEND 

IF  (LCHECK)  WRITE(9,*)J,X(J,1, LSURF) 

END  IF 

CALCULATE  THE  LIMITS  OF  THE  DOMAIN 


L  =  LSURF 

IF  (LCHECK)  WRITE  (9,*)  ’LIMITS  YMAX 
DO  150  J  =  1,NJ 

XSQ  =  X(J,1, LSURF)  *  X(J,1, LSURF) 


ZMIN  ZMAX' 


L 
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XCU  =  XSQ  *  X(J,1,LSURF) 

FAC1  =  3.0  *  XSQ  /  XTOTSQ 
FAC2  =  2.0  *  XCU  /  XTOTCU 

YTEMP  =  (YMAX2  -  YMAX1)  *  (FAC1  -  FAC2)  +  YMAX1 
ZTEMPB  =  (ZMAX2B  -  ZMAX1B)  *  (FAC1  -  FAC2)  +  ZMAX1B 
ZTEMPT  =  (ZMAX2T  -  ZMAX1T)  *  (FAC1  -  FAC2)  +  ZMAX1T 
DO  148  L=1,NL 
Y(J,NK,L)  =  YTEMP 

148  CONTINUE 
DO  149  K  =  1,NK 

Z(  J,K,  1  )=ZTEMPB 
Z(J,K,NL)=ZTEMPT 

149  CONTINUE 

IF  (LCHECK)  WRITE  (9,151)  J,X(J,1,1),Y(J,NK,1), 

&  Z(J,1,1),Z(J,1,NL) 

150  CONTINUE 

151  FORMAT  (I3,4E12.5) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
Y  DIRECTION 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

DISTRIBUTE  POINTS  ALONG  THE  Y  AXIS  AT  X  =  XEND 

IF  (LCHECK)  WRITE  (9  *)  *  Y  DIRECTION’ 

IF(.NOT.LCRANK)THEN 
J  =  JEND 
L  =  LSURF 
K  =  1 

Y(J,K,L)  =  0.0 

IF  (LCHECK)  WRITE  (9,*)  K,Y(J,K,L) 

INNER  SECTION  OF  WING 

YMID  =  BASE  *  .10 
NUMB  =  KMID-1 
DIST  =  YMID 

DETA  =1.0/  FLOAT(NUMB) 

DY  =  0.04  *  DIST 

CALL  FINDC(DY,DIST,NUMB  ,C,  1) 

DO  180  K  =  2.KMID-1 
KK  =  K  -  1 

ETA  =  DETA  *  FLOAT(KK) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  - 1.0) 

Y(J,K,L)  =  FAC  *  DIST 
IF  (LCHECK)  THEN 
Q  =  Y(J,K,L)  /  DY 

PCT  =  (Y(J,K,L)  -  Y(J,K-1,L))  /  BASE 
WRITE  (9,*)  K, Y (J,K,L),Q,PCT 
END  IF 
180  CONTINUE 


non 


c 

K=KMID 
Y(J,K,L)  =  YMID 

IF(LCHECK)  WRITE(9,*)K,Y(J,K,L) 

C 

C  OUTER  SECTION  OF  THE  WING 
C 

NUMB  =  KEDGE  -  KMID 
DIST  =  BASE -YMID 
DETA  =1.0/  FLOAT(NUMB) 

DY  =  1.0  *  DELTA  /  TAN(OMEGA) 

CALL  FINDC(DY,DIST,NUMB  ,C,  1 ) 

DO  190  K=KMJD+1,  KEDGE- 1 
KK  =  KEDGE  -  K 
ETA  =  DETA  *  FLOAT(KK) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1.0) 

Y(J,K,L)  =  BASE  -  FAC  *  DIST 
IF  (LCHECK)  THEN 
Q  =  (BASE  -  Y(J,K,L))  /  DY 
PCT  =  (Y(J,K,L)  -  Y(J,K-1,L))  /  BASE 
WRITE  (9,*)  K,Y(J,K,L),Q,PCT 
END  IF 

190  CONTINUE 
C 

K  =  KEDGE 
Y(J,K,L)  =  BASE 

IF  (LCHECK)  WRITE  (9,*)  K,Y(J,K,L) 

IF  (LSTOP)  CALL  EXIT 
C 

C  BASED  ON  THE  Y  VALUES  AT  X  =  XEND  DISTRIBUTE  Y  VALUES  FROM 
C  X  =  XSTART  TO  XEND  FOR  K  LESS  THAN  OR  EQUAL  TO  KEDGE 
C 

DO  210  J  =  JSTART+UEND 
FAC  =  (X(J,1,LSURF)  -  XSTART)  /  LEN 
DO  210  K  =  2,KEDGE 
YTEMP  =  Y(JEND,K,LSURF)  *  FAC 
DO  210  L  =  1,NL 
Y(J,K,L)  =  YTEMP 
210  CONTINUE 
C 

C  FOR  THE  REGION  BETWEEN  0  AND  XSTART  WITH  K  LESS  THAN  KEDGE, 
SET  ALL 

C  Y  TO  THE  VALUE  AT  JSTART+1 
C 

JP  =  JSTART  +  1 
DO  220  J  =  1JSTART 
DO  220  K  =  2,KEDGE 
DO  220  L  =  1,NL 
Y(J,K,L)  =  Y(JP,K,LSURF) 

220  CONTINUE 
C 

C  DO  THE  CRANKED  WING 
C 
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ELSE 


DISTRIBUTE  POINTS  ALONG  JEND 


J  =  JEND 
L  =  LSURF 
K  =  1 

Y(J,K,L)  =  0.0 

IF  (LCHECK)  WRITE  (9,*)  K,Y(J,K,L) 

INNER  SECTION  OF  WING 

YMID2  =  BASE2  *  0.08 
NUMB  =  KMID-1 
DIST  =  YMID2 
DETA  =1.0/  FLOAT(NUMB) 

DY  =  0.10*  DIST 

CALL  FINDC(D  Y,DIST,NUMB  ,C,  1 ) 

DO  235  K  =  2,KMID-1 
KK  =  K  -  1 

ETA  =  DETA  *  FLOAT(KK) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1.0) 
Y(J,K,L)  =  FAC  *  DIST 
IF  (LCHECK)  THEN 
Q  =  Y(J,K,L)  /  DY 

PCT  =  (Y(J,K,L)  -  Y(J,K-1,L))  /  BASE2 
WRITE  (9,*)  K,Y(J,K,L),Q,PCT 
END  IF 

135  CONTINUE 


K=KMID 

Y(J,K,L)  =  YMID2 
IF(LCHECK)  WRITE(9  *)K,Y(J,K,L) 

OUTER  SECTION  OF  THE  WING 

NUMB  =  KEDGE  KMID 
DIST  =  BASE2  -  YMID2 
DETA  =1.0/  FLOAT(NUMB) 

DY  =  1 .0  *  DELTA  /  TAN(OMEGA) 

CALL  FINDC(DY,DIST,NUMB  ,C,  1 ) 

DO  240  K  =  KMID+1  .KEDGE- 1 
KK  =  KEDGE  -  K 
ETA  =  DETA  *  FLOAT(KK) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  - 1.0) 
Y(J,KX)  =  BASE2  -  FAC  *  DIST 
IF  (LCHECK)  THEN 
Q  =  (BASE2  -  Y(J,K,L))  /  DY 
PCT  =  (Y(J,K,L)  -  Y(J,K-1,L))  /  BASE2 
WRITE  (9,*)  K,Y(J,K,L),Q,PCT 
ENDIF 
240  CONTINUE 
C 
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K  =  KEDGE 
Y  (J,K,L)  =  BASE2 

IF  (LCHECK)  WRITE  (9,*)  K,Y(J,K,L) 

C 

C  BASED  ON  THE  Y  VALUES  AT  X  =  XEND  DISTRIBUTE  Y  VALUES  FROM 
C  X  =  XSTART  TO  XEND  FOR  K  LESS  THAN  OR  EQUAL  TO  KEDGE 
C 

DO  244  J  =  JSTART+1 JEND 
IF(J.LE.JCRANK)  THEN 

Y  (J,KEDGE,LSURF)=(X(J,  1  ,LSURF)-XST  ART)*TAN(S  WEEP1) 

ELSE 

Y(J,KEDGE,LSURF)=(XCRANK-XSTART)*TAN(SWEEP1)+ 

1  (X(J,  1,LSURF)-XCRANK)*TAN(SWEEP2) 

END  IF 

FAC=Y(J,KEDGE,LSURF)/Y(JEND,KEDGE,LSURF) 

DO  244  K  =  2, KEDGE 
YTEMP  =  Y (JEND,K,LSURF)  *  FAC 
DO  244  L  =  1,NL 
Y(J,K,L)  =  YTEMP 
244  CONTINUE 
C 

C  FOR  THE  REGION  BETWEEN  0  AND  XSTART  WITH  K  LESS  THAN  KEDGE, 
SET  ALL 

C  Y  TO  THE  VALUE  AT  JSTART+I 
C 

JP  =  JSTART  +  1 
DO  246  J  =  1  JSTART 
DO  246  K  =  2.KEDGE 
DO  246  L  =  1,NL 
Y(J,K,L)  =  Y(JP,K,LSURF) 

246  CONTINUE 
END  IF 
C 

C  SET  Y  FROM  KEDGE  TO  NK 
C 

DO  250  J  =  1  ,NJ 

DIST  =  Y(J,NK,LSURF)  -  Y (J,KEDGE,LSURF) 

DY  =  Y ( J,KEDGE,LSURF)  -  Y (J,KEDGE- 1  ,LSURF) 

NUMB  =  NK  -  KEDGE 
DETA  =  1.0/  FLOAT(NUMB) 

CALL  FINDC(D  Y  ,DIST,NUMB  ,C,  1 ) 

DO  250  K  =  KEDGE+l.NK 
KK  =  K  -  KEDGE 
ETA  =  DETA  *  FLOAT(KK) 

FAC  =  (EXP(C*ETA)  - 1.0)  /  (EXP(C)  -  1.0) 

YTEMP  =  Y (J,KEDGE,LSURF)  +  FAC  *  DIST 
DO250L=  1,NL 
Y(J,K,L)  =  YTEMP 
250  CONTINUE 
C 

C  DEFINE  MAX  AND  MIN  PLANES  IN  THE  Z  DIRECTION 
C 

DO  260  J=1,NJ 
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DO  260  K=1,NK 

Z(J,K,1)=Z(J,K,1)*((  1  .-(CURVE*  Y(J,K,  1 )/ 

1  Y(J,NK,1))**2)**.5) 

Z(J,K,NL)=Z(  J,K,NL)*((  1  .-(CURVE*  Y(J,K,NL)/ 

1  Y(J,NK,NL))**2)**.5) 

260  CONTINUE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  Z  DIRECTION 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  CALCULATE  THE  Z  VALUES  ABOVE  THE  SURFACE  FROM  THE  ORIGIN 
C  TO  XTOTAL 
C 

NUMB  =  NL  -  LSURF 
DETA  =1.0/  FLOAT(NUMB) 

FAC1  =  DELTA  *  (X(JSTART+ 1,1, LSURF)  -  XSTART)  /  LEN 
DO  280  J=  1,NJ 
IF  (J  .LE.  JSTART)  THEN 
DZ  =  FAC1 

ELSEIF  (J  .LT.  JEND)  THEN 
FAC  =  (X(J  1, LSURF)  -  XSTART)  /  LEN 
DZ  =  DELTA*  FAC 
ELSE 

DZ  =  DELTA 
END  IF 

DO  280  K=1,NK 
DIST =Z(  J,K,NL) 

CALL  FINDC(DZ,DIST,NUMB  ,C,  1 ) 

DO  280  L  =  LSURF+1,NL 
LL  =  L  -  LSURF 
ETA  =  DETA  *  FLOAT(LL) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1.0) 

Z(J,K,L)=  FAC  *  DIST 
280  CONTINUE 
IF  (LCHECK)  THEN 
WRITE  (9,*)  ’  Z  DIRECTION 
J  =  JEND 
K  =  1 

DZ  =  Z(J,K,LSURF+1)  -  Z(J,K,LSURF) 

DO  285  L  =  LSURF, NL 
Q  =  Z(J,K,L)  /  DZ 
WRITE  (9  *)  LZ(J,K,L),Q 
285  CONTINUE 
END  IF 
C 

C  NOW  FOR  THE  BOTTOM 
C 

S2  =  TAN(OMEGA) 

NUMB  =  LSURF  -  2 
DETA  =1.0/  FLO  AT(NUMB) 


nno 


b> 


» 
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DO  290  J  =  JSTART+1 JEND 
DO  290  K  =  1  ,NK 
C 

C  FIRST  LOCATE  THE  Z  VALUE  OF  THE  BOTTOM  SURFACE 
C 

IF(K.LT.KEDGE)THEN 
B2  =  -S2  *  Y (J,KEDGE,LSURF) 

Z(J,K,LSURF-1)  =  Y(J,K,LSURF)  *  S2  +  B2 
ELSE 

Z(J,K,LSURF-1)=Z(J,KEDGE- 1  ,LSURF- 1 ) 

END  IF 
C 

C  CALCULATE  THE  LOCAL  DELTA 
C 

FAC  =  (X(J,1, LSURF)  -  XSTART)  /LEN 
DZ  =  DELTA*  FAC 

DIST  =  ABS(Z(J,K,I)  -  Z(J,K, LSURF- 1)) 

CALL  FINDC(DZ,DIST,NUMB,C,1) 

DO  290  L  =  I,LSURF-2 
LL  =  LSURF  -  L  -  1 
ETA  =  DETA  *  FLOAT(LL) 

FAC  =  (EXP(C*ETA)  -  1.0)  /  (EXP(C)  -  1 .0) 

Z(J,K,L)=Z(  J,K,LS  URF- 1  )-FAC*DIST 
290  CONTINUE 
C 

C  REGION  BETWEEN  0  AND  XSTART,  Z  VALUES  ARE  THOSE  AT  XSTART+1 
C 

JP  =  JSTART  +  1 

Z1  =  Z(JP,1,LSURF+1)  -  Z(JP,1, LSURF) 

Z2  =  Z(JP,1, LSURF)  -  Z(JP,1, LSURF- 1) 
n  7  =  72-71 
DO  310  L  =  1, LSURF- 1 
ZTEMP  =  Z(JP,1,L)  +  DZ 
DO  3 10  J  =  1,  JSTART 
DO  310  K  =  1,NK 
Z(J,K,L)  =  ZTEMP 
310  CONTINUE 
C 

C  OUTPUT  LOOPS 
C 

C  ARC3D  LOOP 
C 

WRITE  (10)  NJ,NK,NL 
DO  999  L  =  1  ,NL 

WRITE  (10)  ((X(J,K,L)J  =  1,NJ),K  =  1,NK), 

*  0Y(J,K,L)J  =  1  ,NJ),K  =  1,NK), 

*  ((Z(J,K,L),J  =  1,NJ),K  =  1,NK) 

999  CONTINUE 

C 

C  WRITE  OUT  GRID  PLANES  FOR  XYGRID  PLOTTING 
C 

OP£N(l  l,FILE=’JPLANE’,STATUS='NEW’,FORM='FORMATTED') 
OPEN(12,FILE='KPLANE',STATUS='NEW',FORM='FORMATTED') 
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0PEN(13,FILE='LPLANE',STATUS='NEW',F0RM='F0RMATTED’) 

C 

;=jend 

DO  1000  K=1,NK 
DO  1000  L=1,NL 

1000  WRITE(11,1030)Y(J,K,L),Z(J,K,L) 

K=1 

DO  1010  J=1,NJ 
DO  1010  L=1,NL 

1010  WRTTE( 12,1 030)X(  J,K,L),Z(J,K,L) 

L=LSURF 
DO  1020  J=1,NJ 
DO  1020  K=1,NK 

1020  WRITE(  13, 1030)X(J,K,L),Y  (J,K,L) 

1030  FORMAT(2E15,7) 

C 

C  IRIS  LOOP 

C  WRITE  (10)  NJ,NK,NL 

C  WRITE  (10)  (((X(J,K,L),J  =  1,NJ),K  =  1,NK),L  =  l.NL), 

C  &  (((Y(J,K,L),J  =  1,NJ),K  =  1,NK),L  =  l.NL), 

C  &  (((Z(J,K,L)J  =  1,NJ),K  =  1,NK),L  =  l.NL) 

C 

CALL  EXIT 
C  STOP 
END 

SUBROUTINE  FINDC 

SUBROUTINE  FINDC(DELTA,DIST,ICELL,C,IDIV) 

C  =  2.00 

ETA  =  FLOAT(IDIV)  /  FLOAT(ICELL) 

DIV  =  IDIV 
DO  1001=  1,20 
EC  =  EXP(C) 

ECM1  =  EC  -  1.0 
ECD  =  EXP(C  *  ETA) 

ECDM1  =  ECD  -  1.0 

BOTTOM  =  ECM1  *  ECD  *  ETA  -  ECDM1  *  EC 
TOP  =  DELTA  -  DIST  *  ECDM1  /  ECM1 
C  =  C  +  TOP  /  BOTTOM  /  DIST  *  ECM1  *  ECMI 
100  CONTINUE 
RETURN 
END 
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